# remove everything that we stored in the environment
rm(list=ls())
# Load necessary libraries
library(readr)
library(dplyr)
library(rstatix)
library(gtsummary)
library(ggplot2)
library(ggpubr)
library(openxlsx)
library(gt)
library(stats)
library(grDevices)
library(car)
library(stringr)
library(lubridate)
library(Hmisc)
library(writexl)
library(tidyr)
library(patchwork)
### 0.0 PREPARATION -----------------------------------------------------------
# load prepared and cleaned file from 00_preparation: df_yyyymmdd.csv
setwd("/Users/nataschastoffel/Documents/GitHub/interoception_NS")
df_RRST <- read.csv("/Users/nataschastoffel/Documents/GitHub/interoception_NS/data/processed/05-2025/Data.RRST_FINAL_20250509.csv")
# factor variables
factor_variables <- c('group', 'sex', 'psychotropic_medication')
# Apply as.numeric to the selected columns
df_RRST[factor_variables] <- lapply(df_RRST[factor_variables], as.factor)
# numeric variables
numeric_variables <- c('age', 'bmi', 'bdi', 'stai_s', 'stai_t', 'maia_total',
'sdq', 'ctq_total', 'ias', 'rrst_sensitivity', 'rrst_slope', 'average_RT', 'metascore',
'unpleasantness', 'asthma', 'dizziness', 'breathlessness',
'correct_RT', 'incorrect_RT')
# Apply as.numeric to the selected columns
df_RRST[numeric_variables] <- lapply(df_RRST[numeric_variables], as.numeric)
df <- df_RRST
df$maia_commonfactor <- rowSums(df[, c("maia_note", "maia_attreg", "maia_aware",
"maia_sfreg", "maia_body", "maia_trust")])RRST : Analysis - sample excluding outlier
0) Data and Preparation
Respiratory Interoception; Analysis of participants of site CH from 15.10.23 until 10.09.2024;
1) Having a Look at the Data
1.1) Identifying the Outlier to exclude
# Testing Outlier using boxplot (standard deviation)
df %>%
group_by(group) %>%
identify_outliers(rrst_sensitivity)# A tibble: 1 × 85
group pcode date age sex handedness diagnosis medication.x smoke_crf
<fct> <int> <chr> <dbl> <fct> <int> <int> <int> <int>
1 1 63 08.05.2024 38 1 1 5 1 0
# ℹ 76 more variables: drugs_crf <int>, bmi <dbl>, bdi <dbl>, stai_t <dbl>,
# stai_s <dbl>, ias <dbl>, maia_note <dbl>, maia_distr <dbl>,
# maia_worry <dbl>, maia_attreg <dbl>, maia_aware <dbl>, maia_sfreg <dbl>,
# maia_body <dbl>, maia_trust <dbl>, maia_total <dbl>, sdq <dbl>,
# contraception <int>, menopause <int>, mc_phase <int>, sfmdrs <int>,
# seiz <int>, cgi <int>, ctq_emoab <int>, ctq_physab <int>, ctq_sexab <int>,
# ctq_emoneg <int>, ctq_physneg <int>, ctq_minim <int>, ctq_total <dbl>, …# one outlier identified = P063 for rrst_sensitivity
#### Box-Plot with Jitter
plotBoxplotGroups <- function(df, var , tit ){
plot <- ggplot(df, aes(x = group, y = var, fill = group)) +
geom_boxplot(show.legend = TRUE) +
labs(y = "Score") +
theme_classic() +
scale_fill_manual(values = c("#868686FF", "#BB4038")) +
scale_x_discrete(labels = c("HC", "FND")) +
guides(fill = "none") +
ggtitle(tit) +
xlab("") +
geom_jitter(shape = 16, position = position_jitter(0.2))
return(plot)
}
p1 <- plotBoxplotGroups(df, (df$rrst_sensitivity), "Respiration Sensitivity") +
stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
p1
# Filter the rows where rrst_sensitivity is greater than 8 and display pcode (visually detected as extreme)
filtered_pcodes <- df[df$rrst_sensitivity > 8, "pcode"]
print(filtered_pcodes) [1] 63# P063
# Calculate mean and standard deviation of rrst_sensitivity
mean_rrst <- mean(df$rrst_sensitivity)
sd_rrst <- sd(df$rrst_sensitivity)
# Define the threshold for filtering (mean + 2.5 * SD)
threshold <- mean_rrst + 2.5 * sd_rrst # 7.07
# Filter the pcodes for rrst_sensitivity greater than the threshold
filtered_pcodes <- df[df$rrst_sensitivity > threshold, "pcode"]
# View the result
filtered_pcodes[1] 63# P063
df_exclusion <- df[df$pcode != "63", ]
p1b <- plotBoxplotGroups(df_exclusion, (df_exclusion$rrst_sensitivity), "Respiration Sensitivity") +
stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
print(p1b)
p1b
# Summary of Change - depending on Exclusion of Outlier
summary(df$rrst_sensitivity) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 2.458 3.280 3.282 4.350 8.752 summary(df_exclusion$rrst_sensitivity) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 2.429 3.277 3.222 4.307 6.242 # for statistical reasons we EXCLUDE P063 from the dataset
df <- df[df$pcode != "63", ]1.2) Testing Distribution of Variables
# Normality of Data -------------------
hist(df$rrst_sensitivity) shapiro.test(df$rrst_sensitivity) # normally distributed
hist(df$rrst_slope) # extremely right skeweddf$rrst_slope_log <- log(df$rrst_slope + 1) # Add 1 to avoid log(0) and log to normalize skewed variables
hist(df$rrst_slope_log) # still skewed but better (for visualization)shapiro.test(df$rrst_slope_log) # NOT normally distributed
shapiro.test(df$rrst_slope_log[df$group == "1"]) # patients within themselves are NOT normally distributed
shapiro.test(df$rrst_slope_log[df$group == "0"]) # controls within themselves are normally distributed
hist(df$metascore) shapiro.test(df$metascore)# NOT normally distributed
shapiro.test(df$metascore[df$group == "1"]) # patients within themselves are NOT normally distributed
shapiro.test(df$metascore[df$group == "0"]) # controls within themselves are normally distributed
df$average_RT <- as.numeric(df$average_RT)
hist(df$average_RT) shapiro.test(df$average_RT) # NOT normally distributed
### other variables
shapiro.test(df$age) # normally distributed
shapiro.test(df$bdi) # NOT normally distributed : p < 0.05
shapiro.test(df$bdi[df$group == "1"]) # patients within themselves are normally distributed
shapiro.test(df$bdi[df$group == "0"]) # controls within themselves are NOT normally distributed
shapiro.test(df$stai_t) # normally distributed
shapiro.test(df$stai_s) # NOT normally distributed
shapiro.test(df$stai_s[df$group == "1"]) # patients within themselves are normally distributed
shapiro.test(df$stai_s[df$group == "0"]) # controls within themselves are NOT normally distributed : p < 0.05
shapiro.test(df$ctq_total) # NOT normally distributed : p < 0.05
shapiro.test(df$ctq_total[df$group == "1"]) # patients within themselves are also NOT normally distributed
shapiro.test(df$ctq_total[df$group == "0"]) # controls within themselves are also NOT normally distributed
shapiro.test(df$sdq) # NOT normally distributed : p < 0.05
shapiro.test(df$sdq[df$group == "1"]) # patients within themselves are NOT normally distributed
shapiro.test(df$sdq[df$group == "0"]) # controls within themselves are NOT normally distributed : p < 0.05
shapiro.test(df$maia_total) # normally distributed
shapiro.test(df$maia_commonfactor) # normally distributed
shapiro.test(df$maia_note) # NOT normally distributed
shapiro.test(df$ias) # NOT normally distributed
shapiro.test(df$bmi) # NOT normally distributed with p < 0.05
shapiro.test(df$bmi[df$group == "1"]) # within patients bmi is normally distributed
shapiro.test(df$bmi[df$group == "0"]) # within control bmi is NOT normally distributed with p < 0.051.3) Testing Collinearity of Variables
## testing co-dependency of variables by creating a correlation matrix -----
corr_matrix <- cor(df[, c("age","bmi", "bdi", "stai_s", "stai_t", "ctq_total",
"ias", "maia_total", "sdq", "rrst_sensitivity",
"rrst_slope", "metascore")], use = "complete.obs")
# Set the cutoff for high correlations, that we then want to exclude
cutoff <- 0.7
# Find the indices of correlations that are greater than 0.7 and exclude self-correlations (diagonal)
high_corr <- which(abs(corr_matrix) > cutoff & abs(corr_matrix) < 1, arr.ind = TRUE)
# Display the pairs of variables with high correlations
high_corr_pairs <- data.frame(
Var1 = rownames(corr_matrix)[high_corr[, 1]],
Var2 = colnames(corr_matrix)[high_corr[, 2]],
Correlation = corr_matrix[high_corr]
)
# View the pairs
print(high_corr_pairs) # stai and bdi are highly correlated, which we do accept though as this is a common clinical happening
cor.test(df$bdi, df$stai_t) #
df <- df %>%
mutate(anx_dep_SUM = stai_t + bdi) # create one variable for the trait anxiety and depression to use as one variable in the model2) Group Differences
### Group Differences as Baseline of Variables of Interst
# parametric t-test for NORMALLY distributed data
res<-t.test(df$rrst_sensitivity[df$group=="1"], df$rrst_sensitivity[df$group=="0"])
res # SIG difference between group
Welch Two Sample t-test
data: df$rrst_sensitivity[df$group == "1"] and df$rrst_sensitivity[df$group == "0"]
t = -2.1801, df = 73.97, p-value = 0.03243
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-1.2343241 -0.0554876
sample estimates:
mean of x mean of y
2.881984 3.526890 res<-t.test(df$stai_t[df$group=="1"], df$stai_t[df$group=="0"])
res # SIG difference between groups with p < 0.05
Welch Two Sample t-test
data: df$stai_t[df$group == "1"] and df$stai_t[df$group == "0"]
t = 3.9112, df = 82.97, p-value = 0.0001871
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
4.563215 14.006552
sample estimates:
mean of x mean of y
46.53488 37.25000 res<-t.test(df$maia_total[df$group=="1"], df$maia_total[df$group=="0"])
res # SIG difference between groups with p < 0.05
Welch Two Sample t-test
data: df$maia_total[df$group == "1"] and df$maia_total[df$group == "0"]
t = -3.679, df = 75.258, p-value = 0.0004374
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-6.547208 -1.947651
sample estimates:
mean of x mean of y
19.87420 24.12163 # non paramatric wilcox test for NOT normally distributed data
res<-wilcox.test(df$metascore[df$group=="1"], df$metascore[df$group=="0"])
res # no sig difference between groups
Wilcoxon rank sum exact test
data: df$metascore[df$group == "1"] and df$metascore[df$group == "0"]
W = 1015, p-value = 0.8962
alternative hypothesis: true location shift is not equal to 0res<-wilcox.test(df$rrst_slope[df$group=="1"], df$rrst_slope[df$group=="0"])
res # no sig difference between groups (p = 0.0629)
Wilcoxon rank sum exact test
data: df$rrst_slope[df$group == "1"] and df$rrst_slope[df$group == "0"]
W = 1267, p-value = 0.06209
alternative hypothesis: true location shift is not equal to 0res<-wilcox.test(df$rrst_slope_log[df$group=="1"], df$rrst_slope_log[df$group=="0"])
res # still same result; confirming the data is the same
Wilcoxon rank sum exact test
data: df$rrst_slope_log[df$group == "1"] and df$rrst_slope_log[df$group == "0"]
W = 1267, p-value = 0.06209
alternative hypothesis: true location shift is not equal to 0res<-wilcox.test(df$ias[df$group=="1"], df$ias[df$group=="0"])
res # SIG difference between groups with p < 0.05
Wilcoxon rank sum test with continuity correction
data: df$ias[df$group == "1"] and df$ias[df$group == "0"]
W = 734.5, p-value = 0.01813
alternative hypothesis: true location shift is not equal to 0res<-wilcox.test(df$maia_note[df$group=="1"], df$maia_note[df$group=="0"])
res # no sig difference between groups
Wilcoxon rank sum test with continuity correction
data: df$maia_note[df$group == "1"] and df$maia_note[df$group == "0"]
W = 986.5, p-value = 0.7191
alternative hypothesis: true location shift is not equal to 0res<-wilcox.test(df$bdi[df$group=="1"], df$bdi[df$group=="0"])
res # SIG difference between groups with p < 0.05
Wilcoxon rank sum test with continuity correction
data: df$bdi[df$group == "1"] and df$bdi[df$group == "0"]
W = 1770, p-value = 4.216e-09
alternative hypothesis: true location shift is not equal to 0res<-wilcox.test(df$stai_s[df$group=="1"], df$stai_s[df$group=="0"])
res # SIG difference between groups with p < 0.05
Wilcoxon rank sum test with continuity correction
data: df$stai_s[df$group == "1"] and df$stai_s[df$group == "0"]
W = 1585.5, p-value = 1.079e-05
alternative hypothesis: true location shift is not equal to 0res<-wilcox.test(df$ctq_total[df$group=="1"], df$ctq_total[df$group=="0"])
res # SIG difference between groups with p < 0.05
Wilcoxon rank sum test with continuity correction
data: df$ctq_total[df$group == "1"] and df$ctq_total[df$group == "0"]
W = 1434.5, p-value = 0.001385
alternative hypothesis: true location shift is not equal to 0res<-wilcox.test(df$sdq[df$group=="1"], df$sdq[df$group=="0"])
res # SIG difference between groups with p < 0.05
Wilcoxon rank sum test with continuity correction
data: df$sdq[df$group == "1"] and df$sdq[df$group == "0"]
W = 1838.5, p-value = 1.379e-10
alternative hypothesis: true location shift is not equal to 0res<-wilcox.test(df$ias[df$group=="1"], df$ias[df$group=="0"])
res # no sig difference between groups
Wilcoxon rank sum test with continuity correction
data: df$ias[df$group == "1"] and df$ias[df$group == "0"]
W = 734.5, p-value = 0.01813
alternative hypothesis: true location shift is not equal to 0res<-wilcox.test(df$bmi[df$group=="1"], df$bmi[df$group=="0"])
res # no sig difference between groups
Wilcoxon rank sum test with continuity correction
data: df$bmi[df$group == "1"] and df$bmi[df$group == "0"]
W = 1326.5, p-value = 0.01943
alternative hypothesis: true location shift is not equal to 0#### get mean/median for the main outcome variables per group
mean(df$rrst_sensitivity[df$group=="1"])[1] 2.881984mean(df$rrst_sensitivity[df$group=="0"])[1] 3.52689sd(df$rrst_sensitivity[df$group=="1"])[1] 1.619366sd(df$rrst_sensitivity[df$group=="0"])[1] 1.128228median(df$rrst_slope[df$group=="1"])[1] 63.2525median(df$rrst_slope[df$group=="0"])[1] 32.4818IQR(df$rrst_slope[df$group=="1"])[1] 2880.422IQR(df$rrst_slope[df$group=="0"])[1] 49.09485median(df$rrst_slope_log[df$group=="1"])[1] 4.162821median(df$rrst_slope_log[df$group=="0"])[1] 3.510839IQR(df$rrst_slope_log[df$group=="1"])[1] 4.884927IQR(df$rrst_slope_log[df$group=="0"])[1] 1.484032median(df$metascore[df$group=="1"])[1] 0.5852489median(df$metascore[df$group=="0"])[1] 0.5727064IQR(df$metascore[df$group=="1"])[1] 0.07714244IQR(df$metascore[df$group=="0"])[1] 0.054550832.1) Demographic Summary Table
# define required subgroups
df_female <- df[df$sex == "1", ]
df_male <- df[df$sex == "2", ]
df_HC <- df[df$group == "0", ]
df_FND <- df[df$group == "1", ]
# overview across group and sex distribution
tbl <- df %>%
mutate(sex = factor(sex, levels = c("1", "2"), labels = c("female", "male"))) %>%
mutate(group = factor(group, levels = c("0", "1"), labels = c("HC", "FND")))
tbl1<-table(tbl$group,tbl$sex)
tbl1
female male
HC 35 13
FND 31 12# overview for mean age across group
summary_stats <- tbl %>%
group_by(group) %>%
get_summary_stats(age, type = "full")
print(summary_stats)# A tibble: 2 × 14
group variable n min max median q1 q3 iqr mad mean sd
<fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 HC age 48 18 67 37 28 50 22 14.1 37.8 13.0
2 FND age 43 18 64 39 31.5 47 15.5 11.9 38.5 11.8
# ℹ 2 more variables: se <dbl>, ci <dbl>res<-t.test(df$age[df$group == "0"], df$age[df$group == "1"]) #
res # no sig difference between the groups in age
Welch Two Sample t-test
data: df$age[df$group == "0"] and df$age[df$group == "1"]
t = -0.25959, df = 88.975, p-value = 0.7958
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-5.849281 4.497537
sample estimates:
mean of x mean of y
37.81250 38.48837 # SUMMARY TABLE DEMOGRAPHICS -----
df_tbl <- df %>%
mutate(sex = factor(sex, levels = c("1", "2"), labels = c("male", "female"))) %>%
mutate(group = factor(group, levels = c("0", "1"), labels = c("HC", "FND"))) %>%
mutate(psychotropic_medication = factor(psychotropic_medication, levels = c("0", "1"), labels = c("no", "yes")))
# Select variable sfor table
Data.t <- select(df_tbl, c("group", "sex", "age", "psychotropic_medication", "smoke_crf", "bmi", "bdi", "stai_t", "ctq_total", "sdq"))
# Create the summary table
tbl1 <- Data.t %>%
tbl_summary(
by = group,
missing = "no",
statistic = list(
age = "{mean} ({sd})",
sex = "{n} ({p})",
psychotropic_medication = "{n} ({p})",
smoke_crf = "{n} ({p})",
bmi = "{median} ({p25}, {p75})",
bdi = "{median} ({p25}, {p75})",
stai_t = "{mean} ({sd})",
ctq_total = "{median} ({p25}, {p75})",
sdq = "{median} ({p25}, {p75})"
),
digits = list(age = 1, sex = 1, psychotropic_medication = 1, smoke = 1, bmi = 1, bdi = 1, stai_t=1, ctq_total=1, sdq=1),
label = list(age = "Age, mean (SD)",
sex = "Sex, count (%)",
psychotropic_medication= "Intake of psychotropic Medication, count (%)",
smoke_crf = "Smoking, count (%)",
bmi = "Body Mass Index kg/m², median (IQR)",
bdi = "Depression using BDI-II, median (IQR), ",
stai_t = "Trait Anxiety using STAI, median (IQR)",
ctq_total = "Childhood Trauma, using CTQ total, median (IQR)",
sdq = "Dissociation using SDQ-20, median (IQR)")
) %>%
modify_header(list(label ~ "**Variable**")) %>%
modify_caption("**Demographics Overview**") %>%
modify_spanning_header(c("stat_1", "stat_2") ~ "**Group**") %>%
add_p( test = list(
age ~ "t.test", # t test for normally distributed variables
bmi ~ "wilcox.test", # wilcos test for not-normally distributed variables
sex ~ "chisq.test", # chi-square test for factor variables
psychotropic_medication ~ "chisq.test", # Fisher's exact test for small samples
smoke_crf ~ "chisq.test", # Fisher's exact test for small samples
bdi ~ "wilcox.test",
stai_t ~ "t.test",
ctq_total = "wilcox.test",
sdq="wilcox.test"
)) %>%
add_overall(last = FALSE)
# Display the table
tbl1| Variable | Overall, N = 911 | Group | p-value2 | |
|---|---|---|---|---|
| HC, N = 481 | FND, N = 431 | |||
| Sex, count (%) | >0.9 | |||
| male | 66.0 (72.5) | 35.0 (72.9) | 31.0 (72.1) | |
| female | 25.0 (27.5) | 13.0 (27.1) | 12.0 (27.9) | |
| Age, mean (SD) | 38.1 (12.4) | 37.8 (13.0) | 38.5 (11.8) | 0.8 |
| Intake of psychotropic Medication, count (%) | 21.0 (23.6) | 2.0 (4.3) | 19.0 (44.2) | <0.001 |
| Smoking, count (%) | 15 (16) | 6 (13) | 9 (21) | 0.4 |
| Body Mass Index kg/m², median (IQR) | 23.7 (21.7, 26.3) | 22.7 (21.7, 25.1) | 25.3 (21.8, 27.8) | 0.019 |
| Depression using BDI-II, median (IQR), | 8.0 (3.0, 15.0) | 4.0 (1.0, 8.0) | 15.0 (9.5, 23.5) | <0.001 |
| Trait Anxiety using STAI, median (IQR) | 41.6 (12.1) | 37.3 (10.3) | 46.5 (12.1) | <0.001 |
| Childhood Trauma, using CTQ total, median (IQR) | 37.0 (30.0, 54.5) | 34.0 (29.0, 41.3) | 46.0 (35.0, 60.5) | 0.001 |
| Dissociation using SDQ-20, median (IQR) | 27.0 (22.0, 36.0) | 23.0 (21.0, 26.0) | 36.0 (28.5, 45.0) | <0.001 |
| 1 n (%); Mean (SD); Median (IQR) | ||||
| 2 Pearson’s Chi-squared test; Welch Two Sample t-test; Wilcoxon rank sum test | ||||
# save table as excel
tbl1_df <- tbl1 %>% as_tibble()
write_xlsx(tbl1_df, "2.1_Table_demographics.xlsx")2.2) FND Summary Table
# SUMMARY TABLE for FND symptoms -----
# create variable of duation of months
df_FND <- df_FND %>%
mutate(duration_months = interval(date_symptom_onset, date) %/% months(1))
as.numeric(df_FND$duration_months) [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
[26] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA# select variables for table
Data.t <- dplyr::select(df_FND, c("sex", "duration_months", "fds", "motor", "weakness", "sensory", "pppd", "cognitive"))
Data.t <- Data.t %>%
mutate(sex = factor(sex, levels = c("1", "2"), labels = c("female", "male")))
# Convert 0/1 values to "no"/"yes" for table, so that it will only show the count "yes"
Data.t <- Data.t %>%
mutate(
across(c(fds, motor, weakness, sensory, pppd, cognitive),
~ factor(.x, levels = c(0, 1), labels = c("no", "yes")))
) %>%
dplyr::select(sex, duration_months, fds, motor, weakness, sensory, pppd, cognitive)
# Create the summary table (which symptoms present with how many patients ; per patient more than one symptom-category possible))
tbl1 <- Data.t %>%
tbl_summary(
by = sex,
missing = "no",
statistic = list(
duration_months = "{mean} ({sd})",
fds = "{n} ({p})",
motor = "{n} ({p})",
weakness = "{n} ({p})",
sensory = "{n} ({p})",
pppd = "{n} ({p})",
cognitive = "{n} ({p})"
),
digits = list(medication = 1, duration_months = 1, fds = 1, motor = 1, weakness = 1, sensory = 1, pppd = 1, cognitive = 1),
label = list(
sex = "Sex",
duration_months = "Symptom Duration [in months]",
fds = "Functional Dissociative Seizures [yes]",
motor = "motor + symptoms [yes]",
weakness = "motor - symptoms [yes]",
sensory = "sensory symptoms [yes]",
pppd = "dizziness (PPPD) [yes]",
cognitive = "cognitive symptoms [yes]"
)
) %>%
modify_header(list(label ~ "**Variable**")) %>%
modify_caption("**Clinical Characteristics**") %>%
modify_spanning_header(c("stat_1", "stat_2") ~ "**Group**") %>%
add_p(
test = list(
duration_months = "t.test",
fds = "fisher.test",
motor = "fisher.test",
weakness = "fisher.test",
sensory = "fisher.test",
ctq_total = "fisher.test",
pppd = "fisher.test",
cognitive = "fisher.test"
)
) %>% add_overall(last = FALSE)
tbl1| Variable | Overall, N = 431 | Group | p-value2 | |
|---|---|---|---|---|
| female, N = 311 | male, N = 121 | |||
| Symptom Duration [in months] | NA (NA) | NA (NA) | NA (NA) | |
| Functional Dissociative Seizures [yes] | 8.0 (18.6) | 8.0 (25.8) | 0.0 (0.0) | 0.082 |
| motor + symptoms [yes] | 20.0 (46.5) | 13.0 (41.9) | 7.0 (58.3) | 0.5 |
| motor - symptoms [yes] | 29.0 (67.4) | 20.0 (64.5) | 9.0 (75.0) | 0.7 |
| sensory symptoms [yes] | 19.0 (44.2) | 15.0 (48.4) | 4.0 (33.3) | 0.5 |
| dizziness (PPPD) [yes] | 2.0 (4.7) | 2.0 (6.5) | 0.0 (0.0) | >0.9 |
| cognitive symptoms [yes] | 4.0 (9.3) | 4.0 (12.9) | 0.0 (0.0) | 0.6 |
| 1 Mean (SD); n (%) | ||||
| 2 Fisher’s exact test | ||||
# save table as excel
tbl1_df <- tbl1 %>% as_tibble()
write_xlsx(tbl1_df, "2.2_Table_clinical_characteristics.xlsx")2.3) Visualization of Group Differences
Here we visualize the different distribution per group across all variables
df$group <- as.factor(df$group)
# plot function for normally dsitributed data, allowing the use of the t test (and the mean)
plotViolinGroups_mean <- function(df, var , tit ){
plot <-ggplot(df, aes(x=group, y=var, fill=group)) +
geom_violin(show.legend = FALSE, width = 0.7, trim = FALSE, alpha = 0.8)+
stat_summary(show.legend = FALSE, fun = mean, geom = "point", color = "black", size = 2) +
labs(y = "Score") +
theme_classic()+
scale_fill_manual(values = c("#868686FF", "#BB4038"))+
scale_x_discrete(labels=c("Controls", "Patients"))+
guides(fill=guide_legend(title="Timepoint"))+
ggtitle(tit) +
xlab("")
return(plot)
}
# plot function for non normally dsitributed data; the wilcox test (and the median)
plotViolinGroups_median <- function(df, var , tit ){
plot <-ggplot(df, aes(x=group, y=var, fill=group)) +
geom_violin(show.legend = FALSE, width = 0.7, trim = FALSE, alpha = 0.8)+
stat_summary(show.legend = FALSE, fun = median, geom = "point", color = "black", size = 2) +
labs(y = "Score") +
theme_classic()+
scale_fill_manual(values = c("#868686FF", "#BB4038"))+
scale_x_discrete(labels=c("Controls", "Patients"))+
ggtitle(tit) +
xlab("")
return(plot)
}
# adding jitter points for visualization in publication
jitteredViolin <- function(df, var, tit){
ggplot(df, aes(x = group, y = {{ var }}, fill = group)) +
geom_violin(show.legend = FALSE, width = 0.7, trim = FALSE, alpha = 0.8) +
geom_jitter(width = 0.1, size = 1.2, color = "black", alpha = 0.6) +
stat_summary(show.legend = FALSE, fun = median, geom = "point", color = "black", size = 2) +
labs(y = "Score") +
theme_classic() +
theme(legend.position = "none") +
scale_fill_manual(values = c("#868686FF", "#BB4038")) +
scale_x_discrete(labels = c("Controls", "Patients")) +
ggtitle(tit) +
xlab("")
}
p1 <- plotViolinGroups_median(df,df$bdi,"BDI - Depression")+ ggpubr::stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
p1
pdf("2.3_Plot_BDI_violin.pdf")
p1
dev.off()quartz_off_screen
2 p2 <- plotViolinGroups_mean(df,df$stai_s,"STAI - State Anxiety")+ stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
p3 <- plotViolinGroups_median(df,df$stai_t,"STAI - Trait Anxiety")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
gridExtra::grid.arrange(p2, p3, ncol = 2)
pdf("2.3_Plot_STAI_violin.pdf")
gridExtra::grid.arrange(p2, p3, ncol = 2)
dev.off()quartz_off_screen
2 p4 <- plotViolinGroups_median(df,df$ctq_total,"CTQ - Childhood Trauma")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
p4
pdf("2.3_Plot_CTQ_violin.pdf")
p4
dev.off()quartz_off_screen
2 p5 <- plotViolinGroups_median(df,df$sdq,"SDQ-20 - Dissociation")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
p5
pdf("2.3_Plot_SDQ-20_violin.pdf")
p5
dev.off()quartz_off_screen
2 p6 <- plotViolinGroups_mean(df,df$maia_total,"MAIA - total score")+ stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
p6common <- plotViolinGroups_mean(df,df$maia_commonfactor,"MAIA - common factor")+ stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
p6noticing <- plotViolinGroups_median(df,df$maia_note,"MAIA subscale - Noticing")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
gridExtra::grid.arrange(p6, p6common, p6noticing, ncol = 3)
pdf("2.3_Plots_MAIA_violin.pdf")
gridExtra::grid.arrange(p6, p6common, p6noticing, ncol = 3)
dev.off()quartz_off_screen
2 p7 <- plotViolinGroups_median(df,df$ias,"IAS - Interoceptive Accuracy")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
p7
pdf("2.3_Plot_IAS_violin.pdf")
p7
dev.off()quartz_off_screen
2 p8 <- plotViolinGroups_mean(df,df$rrst_sensitivity,"Respiratory Sensitivity (RRST alpha)")+ stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
p9 <- plotViolinGroups_median(df,df$rrst_slope_log,"Respiratory Slope (RRST beta)")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
p10 <- plotViolinGroups_median(df,df$metascore,"Respiratory Metacognition")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
gridExtra::grid.arrange(p8, p9, p10, ncol = 2)
pdf("2.3_Plots_RRST_violin.pdf")
gridExtra::grid.arrange(p8, p9, p10, ncol = 2)
dev.off()quartz_off_screen
2 Figure2c <- jitteredViolin(df,df$rrst_sensitivity,"Respiratory Sensitivity (RRST alpha)") +
stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
Figure2c
Figure2d <- jitteredViolin(df,df$rrst_slope_log, "Respiratory Slope (RRST beta)") +
stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
Figure2d
Figure3c <- jitteredViolin(df,df$metascore,"Respiratory Metacognition (AUC)") +
stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
Figure3c
#### save images as svg to merge into a tiff for publication
ggsave("Figure2c.svg",
plot = Figure2c,
dpi = 600,
width = 6,
height = 4,
units = "in",
device = "svg")
ggsave("Figure2d.svg",
plot = Figure2d,
dpi = 600,
width = 6,
height = 4,
units = "in",
device = "svg")
ggsave("Figure3c.svg",
plot = Figure3c,
dpi = 600,
width = 6,
height = 4,
units = "in",
device = "svg")
# Table for dependent Variables: group differences
df_tbl <- df %>%
mutate(group = factor(group, levels = c("0", "1"), labels = c("HC", "FND")))
Data.t <- select(df_tbl, c("group", "maia_total", "maia_commonfactor", "maia_note", "ias", "rrst_sensitivity","rrst_slope", "metascore"))
library(gtsummary)
tbl2 <- Data.t %>%
tbl_summary(
by = group,
missing = "no",
statistic = list(
maia_total = "{mean} ({sd})",
maia_commonfactor = "{mean} ({sd})",
maia_note = "{median} ({IQR})",
ias = "{median} ({IQR})",
rrst_sensitivity = "{mean} ({sd})",
rrst_slope = "{median} ({IQR})",
metascore = "{median} ({IQR})"
),
digits = list(maia_total = 2, maia_commonfactor = 2, maia_note = 2, ias = 2,
rrst_sensitivity = 2, rrst_slope = 2, metascore = 2),
label = list(maia_total= "Interoceptive Awareness Self-Report - MAIA, Mean, (SD)",
maia_commonfactor= "MAIA common factor, Mean, (SD)",
maia_note= "MAIA subscale Noticing, Median, (IQR)",
ias= "Interoceptive Accuracy Self-Report - IAS, Median (IQR)",
rrst_sensitivity= "Respiratory Sensitivity - RRST, Mean (SD)",
rrst_slope= "Interoceptive Decision Prescision - RRST slope, Median (IQR)",
metascore= "Respiratory Metacognition - AUC, Median (IQR)"
)
) %>%
modify_header(list(label ~ "**Variable**")) %>%
modify_caption("**Interoception Overview**") %>%
modify_spanning_header(c("stat_1", "stat_2") ~ "**Group**") %>%
add_p(
test = list(
maia_total = "t.test",
maia_commonfactor = "wilcox.test",
maia_note = "t.test",
ias = "wilcox.test",
rrst_sensitivity = "t.test",
rrst_slope = "wilcox.test",
metascore = "wilcox.test")
) %>%
add_overall(last = FALSE)
tbl2| Variable | Overall, N = 911 | Group | p-value2 | |
|---|---|---|---|---|
| HC, N = 481 | FND, N = 431 | |||
| Interoceptive Awareness Self-Report - MAIA, Mean, (SD) | 22.11 (5.78) | 24.12 (4.49) | 19.87 (6.27) | <0.001 |
| MAIA common factor, Mean, (SD) | 17.15 (5.49) | 19.07 (3.82) | 15.01 (6.27) | <0.001 |
| MAIA subscale Noticing, Median, (IQR) | 3.25 (1.00) | 3.25 (0.75) | 3.25 (1.63) | 0.4 |
| Interoceptive Accuracy Self-Report - IAS, Median (IQR) | 85.00 (18.50) | 89.00 (14.25) | 80.00 (25.50) | 0.018 |
| Respiratory Sensitivity - RRST, Mean (SD) | 3.22 (1.41) | 3.53 (1.13) | 2.88 (1.62) | 0.032 |
| Interoceptive Decision Prescision - RRST slope, Median (IQR) | 39.01 (421.44) | 32.48 (49.09) | 63.25 (2,880.42) | 0.062 |
| Respiratory Metacognition - AUC, Median (IQR) | 0.57 (0.06) | 0.57 (0.05) | 0.59 (0.08) | 0.9 |
| 1 Mean (SD); Median (IQR) | ||||
| 2 Welch Two Sample t-test; Wilcoxon rank sum test; Wilcoxon rank sum exact test | ||||
tbl2_df <- tbl2 %>% as_tibble()
write_xlsx(tbl2_df, "2.3_Table_dependent_variables.xlsx")2.4) Effect Sizes and CI
Here we calculate the effect sizes (cohen’s d or OR) and the confidence intervals of the variables to be reported
# calculation of cohens d for variable of interest
# Define the function
calculate_cohens_d <- function(data, variable, group_var, group1, group2) {
# Subset the data for each group
group1_data <- data[[variable]][data[[group_var]] == group1]
group2_data <- data[[variable]][data[[group_var]] == group2]
# Calculate means and standard deviations
mean1 <- mean(group1_data, na.rm = TRUE)
mean2 <- mean(group2_data, na.rm = TRUE)
sd1 <- sd(group1_data, na.rm = TRUE)
sd2 <- sd(group2_data, na.rm = TRUE)
# Calculate sample sizes
n1 <- sum(!is.na(group1_data))
n2 <- sum(!is.na(group2_data))
# Calculate pooled standard deviation
sd_pooled <- sqrt(((n1 - 1) * sd1^2 + (n2 - 1) * sd2^2) / (n1 + n2 - 2))
# Calculate Cohen's d
cohens_d <- (mean1 - mean2) / sd_pooled
# Return Cohen's d
return(cohens_d)
}
# Calculate Cohen's d for variables of interest
cohens_d_rrst <- calculate_cohens_d(df, "rrst_sensitivity", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_rrst) # cohens d of -0.467 means with an effect size of 0.46 FND has LOWER values than HC[1] -0.4666721cohens_d_slope <- calculate_cohens_d(df, "rrst_slope", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_slope) # cohens d of 0.509 means with an effect size of 0.51 FND has HIGHER values than HC[1] 0.5088956cohens_d_meta <- calculate_cohens_d(df, "metascore", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_meta) # cohens d of -0.128 means with an effect size of 0.128 FND has LOWER values than HC[1] -0.1279049cohens_d_maia <- calculate_cohens_d(df, "maia_total", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_maia) # cohens d of -0.786 means with an effect size of 0.786 FND has LOWER values than HC[1] -0.7864429cohens_d_maia <- calculate_cohens_d(df, "maia_commonfactor", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_maia) # cohens d of -0.793 means with an effect size of 0.793 FND has LOWER values than HC[1] -0.7925178cohens_d_maia_note <- calculate_cohens_d(df, "maia_note", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_maia_note) # cohens d of -0.182 means with an effect size of 0.18 FND has LOWER values than HC[1] -0.1824447cohens_d_ias <- calculate_cohens_d(df, "ias", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_ias) # cohens d of -0.651 means with an effect size of 0.651 FND has LOWER values than HC[1] -0.6510104cohens_d_bdi <- calculate_cohens_d(df, "bdi", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_bdi) # cohens d of 1.382615 means with an effect size of 1.382615 FND has HIGHER values than HC[1] 1.382615cohens_d_stai_t <- calculate_cohens_d(df, "stai_t", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_stai_t) # cohens d of 0.829 means with an effect size of 0.829 FND has HIGHER values than HC[1] 0.8285794cohens_d_stai_s <- calculate_cohens_d(df, "stai_s", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_stai_s) # cohens d of 1.080 means with an effect size of 1.080 FND has HIGHER values than HC[1] 1.080068cohens_d_ctq_total <- calculate_cohens_d(df, "ctq_total", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_ctq_total) # cohens d of 0.662 means with an effect size of 0.662 FND has HIGHER values than HC[1] 0.6621177cohens_d_sdq <- calculate_cohens_d(df, "sdq", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_sdq) # cohens d of 1.500 means with an effect size of 1.500 FND has HIGHER values than HC[1] 1.500201cohens_d_bmi <- calculate_cohens_d(df, "bmi", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_bmi) # cohens d of 0.557 means with an effect size of 0.6621177 FND has HIGHER values than HC[1] 0.5569856# for binary variable we calculate the ODDS RATIO
model <- glm(group ~ psychotropic_medication, data = df, family = binomial)
summary(model)
Call:
glm(formula = group ~ psychotropic_medication, family = binomial,
data = df)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.6061 0.2538 -2.389 0.016912 *
psychotropic_medication1 2.8574 0.7854 3.638 0.000275 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 123.28 on 88 degrees of freedom
Residual deviance: 101.51 on 87 degrees of freedom
(2 observations deleted due to missingness)
AIC: 105.51
Number of Fisher Scoring iterations: 4# Extract CI and Odds Ratio
confint_logit <- confint(model) # CI on the log-odds scale
confint_or <- exp(confint_logit) # Exponentiate to get CI for OR
odds_ratio <- exp(coef(model)) # Calculate Odds Ratios
results <- cbind(odds_ratio, confint_or) # Combine OR and CI
results odds_ratio 2.5 % 97.5 %
(Intercept) 0.5454545 0.3268162 0.8883923
psychotropic_medication1 17.4166659 4.5447406 115.31197542.5) MAIA subscores
maia<-dplyr::select(df, pcode, group, maia_attreg, maia_note, maia_distr, maia_worry, maia_aware, maia_sfreg, maia_body, maia_trust)
# Make long data
maia_long <- gather(maia, key = "subscale", value = "score", -pcode, -group) # make data set long (so that every p code has 6 rows, one for each of the subscores of ctq)
maia_long$subscale<-as.factor(maia_long$subscale)
maia_long$subscale <- factor(maia_long$subscale, levels = rev(levels(maia_long$subscale)))
maia_long$score<-as.numeric(maia_long$score)
maia_long$group<-as.factor(maia_long$group)
maia_long$group <- factor(maia_long$group, levels = c(0, 1), labels = c("HC", "FND"))
# Calculate mean, median, CI, SE
library(dplyr)
alpha=0.05
maia.summary <- maia_long %>%
group_by(group,subscale)%>%
dplyr::summarise(
se = sd(score) / sqrt(length(score)),
t=qt((1-alpha)/2 + .5, length(score)-1),
CI=t*se,
score = mean(score), .groups = 'drop')
maia.summary # according to this order name the lables in the following plot p# A tibble: 16 × 6
group subscale se t CI score
<fct> <fct> <dbl> <dbl> <dbl> <dbl>
1 HC maia_worry 0.144 2.01 0.289 2.85
2 HC maia_trust 0.161 2.01 0.324 3.51
3 HC maia_sfreg 0.151 2.01 0.303 2.82
4 HC maia_note 0.0912 2.01 0.183 3.43
5 HC maia_distr 0.111 2.01 0.223 2.21
6 HC maia_body 0.136 2.01 0.273 2.69
7 HC maia_aware 0.119 2.01 0.240 3.68
8 HC maia_attreg 0.119 2.01 0.239 2.93
9 FND maia_worry 0.166 2.02 0.334 2.97
10 FND maia_trust 0.224 2.02 0.453 2.39
11 FND maia_sfreg 0.201 2.02 0.407 1.88
12 FND maia_note 0.182 2.02 0.367 3.26
13 FND maia_distr 0.178 2.02 0.359 1.90
14 FND maia_body 0.207 2.02 0.417 2.10
15 FND maia_aware 0.214 2.02 0.432 3.08
16 FND maia_attreg 0.176 2.02 0.355 2.30 # STATISTICS
# Multiple t-tests using FDR correction for multiple comparisons
library(rstatix)
stat.test <- maia_long %>%
group_by(subscale) %>%
t_test(score ~ group) %>%
adjust_pvalue(method = "fdr") %>%
add_significance()
stat.test # subscales worrying, trusting, body and attention regulation are sig difference across groups# A tibble: 8 × 11
subscale .y. group1 group2 n1 n2 statistic df p p.adj
<fct> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
1 maia_worry score HC FND 48 43 -0.555 85.7 0.58 0.58
2 maia_trust score HC FND 48 43 4.08 78.0 0.000109 0.000872
3 maia_sfreg score HC FND 48 43 3.73 79.8 0.000353 0.00141
4 maia_note score HC FND 48 43 0.842 62.3 0.403 0.461
5 maia_distr score HC FND 48 43 1.48 71.4 0.144 0.192
6 maia_body score HC FND 48 43 2.40 73.7 0.0189 0.0302
7 maia_aware score HC FND 48 43 2.45 66.4 0.0169 0.0302
8 maia_attreg score HC FND 48 43 2.97 75.0 0.00404 0.0108
# ℹ 1 more variable: p.adj.signif <chr> # Filter for significant results (both * and **)
significant_results <- stat.test %>%
filter(p.adj.signif %in% c("*", "**")) %>%
select(subscale, p.adj.signif) # Keep both subscale and significance level
# VISUALIZATION
# Plot
library(ggsignif)
maia.summary <- maia.summary %>%
left_join(significant_results, by = "subscale") %>%
mutate(significance = ifelse(is.na(p.adj.signif), "", p.adj.signif)) # Add a column for asterisks
dodge <- position_dodge(width=0.9)
maiaplot <- maia.summary %>%
ggplot(aes(y = score, x = subscale, ymin = score - CI, ymax = score + CI, fill = group)) +
geom_bar(position = position_dodge(), stat = "identity") +
geom_errorbar(position = position_dodge(width = 0.9), width = 0.3) +
labs(x = "MAIA Subscales", y = "Scores") +
scale_fill_manual(values = c("#868686FF", "#BB4038")) +
ggtitle("MAIA - Interoceptive Awareness") +
theme_classic(base_size = 9) +
scale_x_discrete(labels = c("Not-Worrying", "Trusting", "Self-Regulation", "Noticing",
"Not-Distracting", "Body Listening", "Emotional Awareness", "Attention Regulation")) +
guides(fill = guide_legend(title = "Groups"))
maiaplot
pdf("2.5_MAIAsubscales.pdf")
maiaplot
dev.off()quartz_off_screen
2 3) Testing potential Confounders for RRST sensitivity
3.2) Anxiety
# RRST Sensitivity GROUP
rrst_affective <- lm(formula = rrst_sensitivity ~ group,
data = df)
summary(rrst_affective) # GROUP
Call:
lm(formula = rrst_sensitivity ~ group, data = df)
Residuals:
Min 1Q Median 3Q Max
-2.88198 -0.94639 -0.07559 0.99151 3.07322
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.5269 0.1995 17.682 <2e-16 ***
group1 -0.6449 0.2902 -2.223 0.0288 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.382 on 89 degrees of freedom
Multiple R-squared: 0.05258, Adjusted R-squared: 0.04194
F-statistic: 4.94 on 1 and 89 DF, p-value: 0.02878confint(rrst_affective, level = 0.95) 2.5 % 97.5 %
(Intercept) 3.130559 3.9232197
group1 -1.221465 -0.0683471# RRST Sensitivity GROUP controlled for affective symptoms
rrst_affective <- lm(formula = rrst_sensitivity ~ group + anx_dep_SUM,
data = df)
summary(rrst_affective) # no sig anymore
Call:
lm(formula = rrst_sensitivity ~ group + anx_dep_SUM, data = df)
Residuals:
Min 1Q Median 3Q Max
-2.94169 -0.91149 -0.05828 1.02851 3.07928
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.715201 0.395404 9.396 6.28e-15 ***
group1 -0.552118 0.336279 -1.642 0.104
anx_dep_SUM -0.004418 0.007999 -0.552 0.582
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.387 on 88 degrees of freedom
Multiple R-squared: 0.05586, Adjusted R-squared: 0.0344
F-statistic: 2.603 on 2 and 88 DF, p-value: 0.07974confint(rrst_affective, level = 0.95) 2.5 % 97.5 %
(Intercept) 2.92941812 4.50098406
group1 -1.22040198 0.11616628
anx_dep_SUM -0.02031377 0.01147804# RRST Sensitivity GROUP controlled for affective symptoms
rrst_anxiety <- lm(formula = rrst_sensitivity ~ group + stai_t,
data = df)
summary(rrst_anxiety) # no sig anymore
Call:
lm(formula = rrst_sensitivity ~ group + stai_t, data = df)
Residuals:
Min 1Q Median 3Q Max
-2.89893 -0.93955 -0.09417 0.97389 3.07027
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.451869 0.529118 6.524 4.21e-09 ***
group1 -0.663605 0.316270 -2.098 0.0388 *
stai_t 0.002014 0.013144 0.153 0.8786
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.39 on 88 degrees of freedom
Multiple R-squared: 0.05284, Adjusted R-squared: 0.03131
F-statistic: 2.454 on 2 and 88 DF, p-value: 0.09177confint(rrst_anxiety, level = 0.95) 2.5 % 97.5 %
(Intercept) 2.40035747 4.50337973
group1 -1.29212600 -0.03508498
stai_t -0.02410788 0.02813585# RRST Sensitivity affective symptoms
rrst_affective <- lm(formula = rrst_sensitivity ~ anx_dep_SUM,
data = df)
summary(rrst_affective) # anx_dep sum score by itself NOT sig for rrst_sensitivity
Call:
lm(formula = rrst_sensitivity ~ anx_dep_SUM, data = df)
Residuals:
Min 1Q Median 3Q Max
-3.2496 -0.8413 0.0435 1.0227 2.8697
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.799079 0.395808 9.598 2.17e-15 ***
anx_dep_SUM -0.010979 0.006995 -1.570 0.12
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.401 on 89 degrees of freedom
Multiple R-squared: 0.02693, Adjusted R-squared: 0.016
F-statistic: 2.463 on 1 and 89 DF, p-value: 0.1201confint(rrst_affective, level = 0.95) 2.5 % 97.5 %
(Intercept) 3.01261753 4.585539931
anx_dep_SUM -0.02487721 0.002919802# RRST Sensitivity anxiety symptoms
rrst_anxiety <- lm(formula = rrst_sensitivity ~ stai_t,
data = df)
summary(rrst_anxiety) # anx_dep sum score by itself NOT sig for rrst_sensitivity
Call:
lm(formula = rrst_sensitivity ~ stai_t, data = df)
Residuals:
Min 1Q Median 3Q Max
-3.2444 -0.7941 0.0841 1.1044 2.9103
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.581435 0.535453 6.689 1.92e-09 ***
stai_t -0.008629 0.012356 -0.698 0.487
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.416 on 89 degrees of freedom
Multiple R-squared: 0.00545, Adjusted R-squared: -0.005725
F-statistic: 0.4877 on 1 and 89 DF, p-value: 0.4868confint(rrst_anxiety, level = 0.95) 2.5 % 97.5 %
(Intercept) 2.51750016 4.64536953
stai_t -0.03317991 0.015922293.2) RRST discomfort
## first we test whether the different control-questions correlate
corr_matrix <- cor(df[, c("breathlessness","dizziness", "unpleasantness", "asthma")], use = "complete.obs")
# Set the cutoff for high correlations, that we then want to exclude
cutoff <- 0.5
# Find the indices of correlations that are greater than 0.7 and exclude self-correlations (diagonal)
high_corr <- which(abs(corr_matrix) > cutoff & abs(corr_matrix) < 1, arr.ind = TRUE)
# Display the pairs of variables with high correlations
high_corr_pairs <- data.frame(
Var1 = rownames(corr_matrix)[high_corr[, 1]],
Var2 = colnames(corr_matrix)[high_corr[, 2]],
Correlation = corr_matrix[high_corr]
)
# View the pairs
print(high_corr_pairs) # no high correlation; no general factor of discomfort[1] Var1 Var2 Correlation
<0 rows> (or 0-length row.names)### calculate correlations to respiratory sensitivity (and control for multiple comparision)
# Run all correlations and store p-values
p_values <- c(
cor.test(df$breathlessness, df$rrst_sensitivity)$p.value,
cor.test(df$asthma, df$rrst_sensitivity)$p.value,
cor.test(df$unpleasantness, df$rrst_sensitivity)$p.value,
cor.test(df$dizziness, df$rrst_sensitivity)$p.value
)
# Apply FDR correction
p_adjusted <- p.adjust(p_values, method = "fdr")
# Optional: associate each result with variable name
results <- data.frame(
variable = c("breathlessness", "asthma", "unpleasantness", "dizziness"),
raw_p = p_values,
p_fdr = p_adjusted)
print(results) variable raw_p p_fdr
1 breathlessness 0.025679985 0.051359970
2 asthma 0.001684408 0.006737632
3 unpleasantness 0.302098918 0.402798557
4 dizziness 0.921408315 0.921408315# Calculate Cohen's d for variables of interest
cohens_d_discomfort <- calculate_cohens_d(df, "breathlessness", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_discomfort) # cohens d of -0.467 means with an effect size of 0.46 FND has LOWER values than HC[1] 0.6465991cohens_d_discomfort <- calculate_cohens_d(df, "asthma", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_discomfort) # cohens d of -0.467 means with an effect size of 0.46 FND has LOWER values than HC[1] 0.4498372cohens_d_discomfort <- calculate_cohens_d(df, "unpleasantness", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_discomfort) # cohens d of -0.467 means with an effect size of 0.46 FND has LOWER values than HC[1] 0.4801304cohens_d_discomfort <- calculate_cohens_d(df, "dizziness", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_discomfort) # cohens d of -0.467 means with an effect size of 0.46 FND has LOWER values than HC[1] 0.4457795## testing linear models for the effect on the identified group difference
# RRST Sensitivity GROUP controlled for discomfort of the task
rrst_discomfort <- lm(formula = rrst_sensitivity ~ group + unpleasantness,
data = df)
summary(rrst_discomfort) # not sig anymore; p = 0.0691
Call:
lm(formula = rrst_sensitivity ~ group + unpleasantness, data = df)
Residuals:
Min 1Q Median 3Q Max
-3.01017 -0.92969 -0.01731 1.06434 2.95843
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.665167 0.271946 13.478 <2e-16 ***
group1 -0.558640 0.303399 -1.841 0.0691 .
unpleasantness -0.003455 0.005862 -0.589 0.5572
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.38 on 85 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.05026, Adjusted R-squared: 0.02791
F-statistic: 2.249 on 2 and 85 DF, p-value: 0.1117confint(rrst_discomfort, level = 0.95) 2.5 % 97.5 %
(Intercept) 3.12446522 4.205868773
group1 -1.16187795 0.044597149
unpleasantness -0.01510882 0.008199755# RRST Sensitivity GROUP controlled for discomfort of the task
rrst_asthma <- lm(formula = rrst_sensitivity ~ group + asthma,
data = df)
summary(rrst_asthma) # group not sig, but asthma with p. = 0.005
Call:
lm(formula = rrst_sensitivity ~ group + asthma, data = df)
Residuals:
Min 1Q Median 3Q Max
-3.2401 -0.6751 0.0155 0.9793 2.7620
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.657422 0.195750 18.684 < 2e-16 ***
group1 -0.417360 0.289355 -1.442 0.15287
asthma -0.022484 0.007853 -2.863 0.00528 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.32 on 85 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.1303, Adjusted R-squared: 0.1098
F-statistic: 6.365 on 2 and 85 DF, p-value: 0.002656confint(rrst_asthma, level = 0.95) 2.5 % 97.5 %
(Intercept) 3.26821773 4.046625986
group1 -0.99267486 0.157954761
asthma -0.03809798 -0.006869975# RRST Sensitivity GROUP controlled for discomfort of the task
rrst_breath <- lm(formula = rrst_sensitivity ~ group + breathlessness,
data = df)
summary(rrst_breath) # both factors not sig
Call:
lm(formula = rrst_sensitivity ~ group + breathlessness, data = df)
Residuals:
Min 1Q Median 3Q Max
-3.1858 -0.7460 -0.0539 0.9959 2.9155
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.701317 0.215022 17.214 <2e-16 ***
group1 -0.437026 0.305481 -1.431 0.1562
breathlessness -0.010756 0.006226 -1.728 0.0877 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.359 on 85 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.07873, Adjusted R-squared: 0.05706
F-statistic: 3.632 on 2 and 85 DF, p-value: 0.03065confint(rrst_breath, level = 0.95) 2.5 % 97.5 %
(Intercept) 3.27379618 4.128838668
group1 -1.04440446 0.170352839
breathlessness -0.02313453 0.001621564# RRST Sensitivity GROUP controlled for discomfort of the task
rrst_dizzy <- lm(formula = rrst_sensitivity ~ group + dizziness,
data = df)
summary(rrst_dizzy) # group still sig p = 0.042
Call:
lm(formula = rrst_sensitivity ~ group + dizziness, data = df)
Residuals:
Min 1Q Median 3Q Max
-2.99484 -0.79845 -0.04541 0.97905 3.02539
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.507924 0.245147 14.310 <2e-16 ***
group1 -0.624278 0.302629 -2.063 0.0422 *
dizziness 0.001906 0.005379 0.354 0.7239
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.382 on 85 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.04778, Adjusted R-squared: 0.02538
F-statistic: 2.133 on 2 and 85 DF, p-value: 0.1248confint(rrst_dizzy, level = 0.95) 2.5 % 97.5 %
(Intercept) 3.020507382 3.99534123
group1 -1.225986511 -0.02257024
dizziness -0.008789153 0.01260153### Visualization
shapiro.test(df$breathlessness)# NOT normally distributed
Shapiro-Wilk normality test
data: df$breathlessness
W = 0.80911, p-value = 2.605e-09shapiro.test(df$dizziness)# NOT normally distributed
Shapiro-Wilk normality test
data: df$dizziness
W = 0.87343, p-value = 4.003e-07shapiro.test(df$asthma)# NOT normally distributed
Shapiro-Wilk normality test
data: df$asthma
W = 0.51855, p-value = 1.662e-15shapiro.test(df$unpleasantness)# NOT normally distributed
Shapiro-Wilk normality test
data: df$unpleasantness
W = 0.95164, p-value = 0.002433discomfort1 <- jitteredViolin(df, df$unpleasantness, "A Unpleasantness") +
stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1, label.y = 110) +
ylim(0, 120)
discomfort2 <- jitteredViolin(df, df$dizziness, "B Dizziness") +
stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1, label.y = 110) +
ylim(0, 120)
discomfort3 <- jitteredViolin(df, df$breathlessness, "C Breathlessness") +
stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1, label.y = 110) +
ylim(0, 120)
discomfort4 <- jitteredViolin(df, df$asthma, "D Asthma") +
stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1, label.y = 110) +
ylim(0, 120)
gridExtra::grid.arrange(discomfort1, discomfort2, discomfort3, discomfort4, ncol = 2, nrow = 2)
pdf("3.2_Plots_discomfort_violin.pdf")
gridExtra::grid.arrange(discomfort1, discomfort2, discomfort3, discomfort4, ncol = 2)
dev.off()quartz_off_screen
2 ggsave("Figure6.tiff",
plot = gridExtra::grid.arrange(discomfort1, discomfort2, discomfort3, discomfort4, ncol = 2),
dpi = 600,
width = 6,
height = 4,
units = "in",
device = "tiff",
compression = "lzw")
3.3) Response-Time
Note: for analysis on Response Time P029 had to be excluded as this patient was not able to press the response button themselves, but rather gave a sign for the experimenter to click with a mouseclick; thus the response time do not mirror properly the response time
# exclude particioant for all analysis on Resposne time and create new df
df_RT <- df[df$pcode != "29", ]
df_RT_HC <- df_RT[df_RT$group == "0", ]
df_RT_FND <- df_RT[df_RT$group == "1", ]
# mean age
df_RT %>%
group_by(group) %>%
get_summary_stats(average_RT, type = "mean_sd") # A tibble: 2 × 5
group variable n mean sd
<fct> <fct> <dbl> <dbl> <dbl>
1 0 average_RT 48 1.14 0.392
2 1 average_RT 42 1.07 0.414median(df$average_RT[df$group=="1"])[1] 0.98547median(df$average_RT[df$group=="0"])[1] 1.109339IQR(df$average_RT[df$group=="1"])[1] 0.5055542IQR(df$average_RT[df$group=="0"])[1] 0.5107342res<-wilcox.test(df$average_RT[df$group=="1"], df$average_RT[df$group=="0"])
res # no sig difference between groups
Wilcoxon rank sum exact test
data: df$average_RT[df$group == "1"] and df$average_RT[df$group == "0"]
W = 912, p-value = 0.3438
alternative hypothesis: true location shift is not equal to 0## first we test whether response time correlates with the sensitivity score (as in, whether or not the sensitivity score is mirroring the impulsiveness or the jumping to conclusion)
cor.test(df_RT$average_RT, df_RT$rrst_sensitivity) # r = -0.26, p = 0.013
Pearson's product-moment correlation
data: df_RT$average_RT and df_RT$rrst_sensitivity
t = -2.5445, df = 88, p-value = 0.01269
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.44476592 -0.05783141
sample estimates:
cor
-0.2617873 cor.test(df_RT$average_RT, df_RT$rrst_slope) # also testing the slope to represent interoceptive precision; p = 0.458
Pearson's product-moment correlation
data: df_RT$average_RT and df_RT$rrst_slope
t = 0.74444, df = 88, p-value = 0.4586
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1301134 0.2815872
sample estimates:
cor
0.07910932 # separately per group Pearson
cor.test(df_RT_FND$average_RT, df_RT_FND$rrst_sensitivity) # p = 0.060
Pearson's product-moment correlation
data: df_RT_FND$average_RT and df_RT_FND$rrst_sensitivity
t = -1.9321, df = 40, p-value = 0.06045
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.54748026 0.01291418
sample estimates:
cor
-0.2921639 cor.test(df_RT_FND$average_RT, df_RT_FND$rrst_slope) # p = 0.3745
Pearson's product-moment correlation
data: df_RT_FND$average_RT and df_RT_FND$rrst_slope
t = 0.99143, df = 40, p-value = 0.3274
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1564263 0.4381747
sample estimates:
cor
0.1548679 cor.test(df_RT_HC$average_RT, df_RT_HC$rrst_sensitivity) # p = 0.041
Pearson's product-moment correlation
data: df_RT_HC$average_RT and df_RT_HC$rrst_sensitivity
t = -2.1061, df = 46, p-value = 0.04068
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.53556671 -0.01356923
sample estimates:
cor
-0.2965605 cor.test(df_RT_HC$average_RT, df_RT_HC$rrst_slope) # p = 0.885
Pearson's product-moment correlation
data: df_RT_HC$average_RT and df_RT_HC$rrst_slope
t = 0.14545, df = 46, p-value = 0.885
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.2643045 0.3037247
sample estimates:
cor
0.02144032 # RRST Sensitivity GROUP controlled for average RT
rrst_rt <- lm(formula = rrst_sensitivity ~ group + average_RT,
data = df_RT)
summary(rrst_rt) # both sig
Call:
lm(formula = rrst_sensitivity ~ group + average_RT, data = df_RT)
Residuals:
Min 1Q Median 3Q Max
-3.2230 -0.9557 -0.0913 1.0190 3.4978
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.6692 0.4472 10.441 < 2e-16 ***
group1 -0.7023 0.2834 -2.478 0.01514 *
average_RT -1.0012 0.3536 -2.831 0.00576 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.336 on 87 degrees of freedom
Multiple R-squared: 0.1299, Adjusted R-squared: 0.1099
F-statistic: 6.497 on 2 and 87 DF, p-value: 0.002345confint(rrst_rt, level = 0.95) 2.5 % 97.5 %
(Intercept) 3.780383 5.5580952
group1 -1.265640 -0.1390279
average_RT -1.704135 -0.2983063# Viszalization
df_plot <- df_RT %>%
mutate(group = factor(group, levels = c("0", "1"), labels = c("HC", "FND")))
RRST_RT<- ggplot(df_plot, aes(x = average_RT, y = rrst_sensitivity, color = as.factor(group))) +
geom_point() +
geom_smooth(method = "lm", se = TRUE) +
labs(x = "Average Response Time in s",
y = "Respiratory Sensitivity (RRST)",
title = "B) Respiratory Interoception and Response Time",
color = "Group") +
theme_classic()+
scale_color_manual(values = c("FND" = "#E95248", "HC" = "#868686FF")) +
theme(panel.grid = element_blank())
RRST_RT
ggsave("Figure5b.tiff",
plot = RRST_RT,
dpi = 600,
width = 6,
height = 4,
units = "in",
device = "tiff",
compression = "lzw")
RTplot <- plotViolinGroups_median(df_RT,df_RT$average_RT,"Average Response Time")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
RTplot
ggsave("Figure5a.tiff",
plot = RTplot,
dpi = 600,
width = 6,
height = 4,
units = "in",
device = "tiff",
compression = "lzw")
######
# Convert to long format
Data.RT <- select(df_RT, c("group", "incorrect_RT", "correct_RT", "pcode")) %>%
mutate(group = factor(group, levels = c("0", "1"), labels = c("HC", "FND")))
df_long <- Data.RT %>%
pivot_longer(cols = c(correct_RT, incorrect_RT),
names_to = "accuracy",
values_to = "RT") %>%
mutate(
accuracy = factor(accuracy, levels = c("correct_RT", "incorrect_RT"),
labels = c("Correct", "Incorrect")),
group = factor(group, levels = c("FND", "HC")),
group_accuracy = paste(group, accuracy, sep = "_")
)
# Define color palette
group_colors <- c("FND" = "#E95248", "HC" = "#868686FF")
fill_colors <- c("FND_Correct" = "#f4cccc",
"FND_Incorrect" = "#E95248",
"HC_Correct" = "#bcbcbc",
"HC_Incorrect" = "#868686FF")
# Within-group comparisons (paired)
within_group_tests <- df_long %>%
group_by(group) %>%
wilcox_test(RT ~ accuracy, paired = TRUE) %>%
add_significance() %>%
mutate(group_accuracy = paste(group, "Correct_vs_Incorrect", sep = "_")) %>%
add_xy_position(x = "group_accuracy", dodge = 0.6)
# Between-group comparisons (unpaired)
between_group_tests <- df_long %>%
group_by(accuracy) %>%
wilcox_test(RT ~ group, paired = FALSE) %>%
add_significance() %>%
mutate(group_accuracy = paste("FND", accuracy, sep = "_")) %>% # use FND_* for x-position
add_xy_position(x = "group_accuracy", dodge = 0.6)
violin_plots_RTs <- ggviolin(df_long, x = "group_accuracy", y = "RT",
fill = "group_accuracy", palette = fill_colors,
add = "jitter", width = 1, trim = FALSE) +
labs(title = "A) Reaction Times by Group and Accuracy",
x = NULL, y = "Reaction Time (RT)") +
theme_classic() +
theme(legend.position = "none")
violin_plots_RTs
ggsave("3.3_violin_plots_RTs.tiff",
plot = violin_plots_RTs,
dpi = 600,
width = 6,
height = 4,
units = "in",
device = "tiff",
compression = "lzw")
ggsave("Figure5.tiff",
plot = gridExtra::grid.arrange(violin_plots_RTs, RRST_RT, ncol = 2),
dpi = 600,
width = 8,
height = 4,
units = "in",
device = "tiff",
compression = "lzw")
ggsave("Figure5a-b.svg",
plot = gridExtra::grid.arrange(violin_plots_RTs, RRST_RT, ncol = 2),
dpi = 600,
width = 8,
height = 4,
units = "in",
device = "svg")
4) Symptom Severity as dependent variables
# Compute the correlation matrix with variables of symptom severity to see whether they correlate
# Select the desired variables
severity_variables <- df_FND[, c("sfmdrs", "cgi", "sdq")]
# Compute the correlation matrix with rcorr
correlation_results <- Hmisc::rcorr(as.matrix(severity_variables))
# Extract correlation coefficients
cor_matrix <- correlation_results$r
# Extract p-values
p_matrix <- correlation_results$P
# Print results
cor_matrix sfmdrs cgi sdq
sfmdrs 1.00000000 0.5912098 0.05993439
cgi 0.59120983 1.0000000 0.30738677
sdq 0.05993439 0.3073868 1.00000000p_matrix sfmdrs cgi sdq
sfmdrs NA 2.986506e-05 0.70262560
cgi 2.986506e-05 NA 0.04495112
sdq 7.026256e-01 4.495112e-02 NA# SIG correlations of different symptom severity scores
### CGI & SFMDRS corr: r= 0.5841928 , p=3.136863e-05
### CGI & SDQ corr: r = 0.31 , p = 0.0454.1 ) FND only
# Exclude participants with sfmdrs = 0
filtered_df <- df_FND[df_FND$sfmdrs != 0, ]
# Count how many participants were excluded
excluded_count <- sum(df_FND$sfmdrs == 0)
excluded_count [1] 10# SIMPLE CORRELATION OF RRST sensitivity and symptom severity scores - including SDQ
cor1 <- cor.test(df_FND$rrst_sensitivity, df_FND$sdq, method = "pearson") # r = -0.3828185, p-value = 0.01129
cor2 <- cor.test(filtered_df$rrst_sensitivity, filtered_df$sfmdrs, method = "pearson")
cor3 <- cor.test(df_FND$rrst_sensitivity, df_FND$cgi, method = "pearson")
p_values <- c(cor1$p.value, cor2$p.value, cor3$p.value)
p_adjusted <- p.adjust(p_values, method = "bonferroni")
results_rrst <- data.frame(
Comparison = c("FND: rrst_sensitivity vs sdq", "FND: rrst_sensitivity vs sfmdrs", "FND: rrst_sensitivity vs cgi"),
Original_P = p_values,
Adjusted_P = p_adjusted)
print(results_rrst) Comparison Original_P Adjusted_P
1 FND: rrst_sensitivity vs sdq 0.01128738 0.03386214
2 FND: rrst_sensitivity vs sfmdrs 0.36151691 1.00000000
3 FND: rrst_sensitivity vs cgi 0.30365832 0.91097497# Save results as a text file
write.table(results_rrst, "3.1_corr_RRST_symptoms.txt", sep = "\t", row.names = FALSE, quote = FALSE)
# metascore
cor1 <- cor.test(df_FND$metascore, df_FND$sdq, method = "pearson") # r = -0.36099, p-value = 0.0174
cor2 <- cor.test(filtered_df$metascore, filtered_df$sfmdrs, method = "pearson")
cor3 <- cor.test(df_FND$metascore, df_FND$cgi, method = "pearson")
p_values <- c(cor1$p.value, cor2$p.value, cor3$p.value)
p_adjusted <- p.adjust(p_values, method = "bonferroni")
results_meta <- data.frame(
Comparison = c("FND: metascore vs sdq", "FND: metascore vs sfmdrs", "FND: metascore vs cgi"),
Original_P = p_values,
Adjusted_P = p_adjusted)
print(results_meta) Comparison Original_P Adjusted_P
1 FND: metascore vs sdq 0.01739608 0.05218824
2 FND: metascore vs sfmdrs 0.62780713 1.00000000
3 FND: metascore vs cgi 0.31527501 0.94582503 # Save results as a text file
write.table(results_rrst, "3.1_corr_meta_symptoms.txt", sep = "\t", row.names = FALSE, quote = FALSE)
#### CLINICAL VARIABLES and INTEROCEPTIVE SENSIBILITY
# separate per group
cor.test(df_FND$maia_total, df_FND$sdq, method = "pearson")
Pearson's product-moment correlation
data: df_FND$maia_total and df_FND$sdq
t = -1.0818, df = 41, p-value = 0.2856
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.4446866 0.1407959
sample estimates:
cor
-0.1665927 cor.test(filtered_df$maia_total, filtered_df$sfmdrs, method = "pearson")
Pearson's product-moment correlation
data: filtered_df$maia_total and filtered_df$sfmdrs
t = 1.1102, df = 31, p-value = 0.2754
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1583914 0.5049599
sample estimates:
cor
0.1955512 cor.test(df_FND$maia_total, df_FND$cgi, method = "pearson")
Pearson's product-moment correlation
data: df_FND$maia_total and df_FND$cgi
t = 0.83381, df = 41, p-value = 0.4092
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1781232 0.4134383
sample estimates:
cor
0.1291288 cor.test(df_FND$ias, df_FND$sdq, method = "pearson") # t = -2.32, df = 41, p-value = 0.02539
Pearson's product-moment correlation
data: df_FND$ias and df_FND$sdq
t = -2.32, df = 41, p-value = 0.02539
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.58150466 -0.04490785
sample estimates:
cor
-0.340657 cor.test(filtered_df$ias, filtered_df$sfmdrs, method = "pearson")
Pearson's product-moment correlation
data: filtered_df$ias and filtered_df$sfmdrs
t = 0.6599, df = 31, p-value = 0.5142
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.2351120 0.4431019
sample estimates:
cor
0.1176971 cor.test(df_FND$ias, df_FND$cgi, method = "pearson")
Pearson's product-moment correlation
data: df_FND$ias and df_FND$cgi
t = -1.602, df = 41, p-value = 0.1168
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.50615656 0.06216586
sample estimates:
cor
-0.2427097 # Collect all p-values from test of interoception varibales and clincial variables
p_values <- c(
cor.test(df_FND$maia_total, df_FND$sdq, method = "pearson")$p.value,
cor.test(filtered_df$maia_total, filtered_df$sfmdrs, method = "pearson")$p.value,
cor.test(df_FND$maia_total, df_FND$cgi, method = "pearson")$p.value,
cor.test(df_FND$ias, df_FND$sdq, method = "pearson")$p.value,
cor.test(filtered_df$ias, filtered_df$sfmdrs, method = "pearson")$p.value,
cor.test(df_FND$ias, df_FND$cgi, method = "pearson")$p.value
)
# Apply Bonferroni correction
adjusted_p_values <- p.adjust(p_values, method = "bonferroni")
# Combine with test names for clarity
test_names <- c(
"MAIA vs sdq (FND)",
"MAIA vs sfmdrs (FND)",
"MAIAvs cgi (FND)",
"IAS vs sdq (FND)",
"IAS vs sfmdrs (FND)",
"IAS vs cgi (FND)")
# Create a data frame for easier interpretation
results <- data.frame(
Test = test_names,
P_Value = p_values,
Adjusted_P_Value = adjusted_p_values)
# View results
print(results) Test P_Value Adjusted_P_Value
1 MAIA vs sdq (FND) 0.2856496 1.0000000
2 MAIA vs sfmdrs (FND) 0.2754433 1.0000000
3 MAIAvs cgi (FND) 0.4092203 1.0000000
4 IAS vs sdq (FND) 0.0253949 0.1523694
5 IAS vs sfmdrs (FND) 0.5141939 1.0000000
6 IAS vs cgi (FND) 0.1168329 0.7009974 # Save results as a text file
write.table(results_rrst, "3.1_corr_intero_symptoms.txt", sep = "\t", row.names = FALSE, quote = FALSE)4.2 ) SDQ and RRST: separate corr per group
cor.test(df_FND$metascore, df_FND$sdq, method = "pearson")
Pearson's product-moment correlation
data: df_FND$metascore and df_FND$sdq
t = -2.4786, df = 41, p-value = 0.0174
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.5966450 -0.0680211
sample estimates:
cor
-0.36099 cor.test(df_FND$rrst_sensitivity, df_FND$sdq, method = "pearson")
Pearson's product-moment correlation
data: df_FND$rrst_sensitivity and df_FND$sdq
t = -2.6534, df = 41, p-value = 0.01129
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.61271421 -0.09318929
sample estimates:
cor
-0.3828185 cor.test(df_HC$metascore, df_HC$sdq, method = "pearson")
Pearson's product-moment correlation
data: df_HC$metascore and df_HC$sdq
t = 0.63293, df = 46, p-value = 0.5299
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1964034 0.3673528
sample estimates:
cor
0.09291647 cor.test(df_HC$rrst_sensitivity, df_HC$sdq, method = "pearson")
Pearson's product-moment correlation
data: df_HC$rrst_sensitivity and df_HC$sdq
t = -0.27916, df = 46, p-value = 0.7814
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.3215025 0.2458833
sample estimates:
cor
-0.04112447 # Prepare the data for visualization
df_plot <- df %>%
mutate(group = factor(group, levels = c(0, 1), labels = c("HC", "FND")))
# Define color palette
group_colors <- c("FND" = "#E95248", "HC" = "#868686FF")
# Panel A
RRST_SDQ <- ggplot(df_plot, aes(x = sdq, y = rrst_sensitivity, color = group)) +
geom_point() +
geom_smooth(method = "lm", se = TRUE) +
labs(
x = "Dissociative Score (SDQ-20)",
y = "Respiratory Sensitivity (RRST)",
title = "A Respiratory Interoception (FND vs HC)",
color = "Group"
) +
theme_classic() +
scale_color_manual(values = group_colors) +
theme(panel.grid = element_blank())
# Panel B
META_SDQ <- ggplot(df_plot, aes(x = sdq, y = metascore, color = group)) +
geom_point() +
geom_smooth(method = "lm", se = TRUE) +
labs(
x = "Dissociative Score (SDQ-20)",
y = "Respiratory Metacognition (AUC)",
title = "B Metacognitive Insight (FND vs HC)",
color = "Group"
) +
theme_classic() +
scale_color_manual(values = group_colors) +
theme(panel.grid = element_blank())
# Combine with patchwork
combined_plot <- RRST_SDQ | META_SDQ # side-by-side layout
# Save as TIFF with 600 DPI
ggsave("Figure4_combined.tiff", plot = combined_plot,
dpi = 600, width = 10, height = 5, units = "in", device = "tiff")`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'4.3 ) Identifying Covariates for Regression Analysis
We identified SDQ to be correlated with RRST sensitivity and metacognition in the FND population only.
Next, we need to identify which covariates will be included in our linear regression (separate per group) to test whether these correlation survive also the control for these covariates. To do so, we choose those variables that are both theoretically but also statistically correlated with both the predictor (RRST sensitivity and metacognition in our case) and the outcome variable (SDQ in our case).
# does age need to be a covariate in the model ?
cor.test(df_FND$age, df_FND$rrst_sensitivity, method = "pearson")
Pearson's product-moment correlation
data: df_FND$age and df_FND$rrst_sensitivity
t = -0.99609, df = 41, p-value = 0.325
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.4340205 0.1537269
sample estimates:
cor
-0.153714 cor.test(df_FND$age, df_FND$metascore, method = "pearson")
Pearson's product-moment correlation
data: df_FND$age and df_FND$metascore
t = 1.0698, df = 41, p-value = 0.291
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1426087 0.4432013
sample estimates:
cor
0.1647936 cor.test(df_FND$age, df_FND$sdq, method = "pearson")
Pearson's product-moment correlation
data: df_FND$age and df_FND$sdq
t = 1.6267, df = 41, p-value = 0.1115
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.05844549 0.50892864
sample estimates:
cor
0.2462205 # does bmi need to be a covariate in the model?
cor.test(df_FND$bmi, df_FND$rrst_sensitivity, method = "pearson")
Pearson's product-moment correlation
data: df_FND$bmi and df_FND$rrst_sensitivity
t = -0.5235, df = 41, p-value = 0.6034
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.3727075 0.2243494
sample estimates:
cor
-0.08148511 cor.test(df_FND$bmi, df_FND$metascore, method = "pearson")
Pearson's product-moment correlation
data: df_FND$bmi and df_FND$metascore
t = -0.58658, df = 41, p-value = 0.5607
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.3811282 0.2150079
sample estimates:
cor
-0.09122709 cor.test(df_FND$bmi, df_FND$sdq, method = "pearson")
Pearson's product-moment correlation
data: df_FND$bmi and df_FND$sdq
t = -0.97486, df = 41, p-value = 0.3353
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.4313571 0.1569247
sample estimates:
cor
-0.1505131 # does anx_dep_sum need to be a covariate in the model?
cor.test(df_FND$anx_dep_SUM, df_FND$rrst_sensitivity, method = "pearson")
Pearson's product-moment correlation
data: df_FND$anx_dep_SUM and df_FND$rrst_sensitivity
t = -1.0783, df = 41, p-value = 0.2872
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.4442492 0.1413301
sample estimates:
cor
-0.1660627 cor.test(df_FND$anx_dep_SUM, df_FND$metascore, method = "pearson")
Pearson's product-moment correlation
data: df_FND$anx_dep_SUM and df_FND$metascore
t = -0.64011, df = 41, p-value = 0.5257
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.3882177 0.2070577
sample estimates:
cor
-0.0994722 cor.test(df_FND$anx_dep_SUM, df_FND$sdq, method = "pearson") # trend with r= 0.294215, p-value = 0.05548
Pearson's product-moment correlation
data: df_FND$anx_dep_SUM and df_FND$sdq
t = 1.9711, df = 41, p-value = 0.05548
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.006722943 0.546285935
sample estimates:
cor
0.294215 # does medication need to be a covariate in the model?
res<-t.test(df_FND$rrst_sensitivity[df_FND$psychotropic_medication == "1"], df_FND$rrst_sensitivity[df_FND$psychotropic_medication=="0"])
res # sig difference across groups with or without intake of psychotropic medicatioN: t = -2.1273, p-value = 0.04011
Welch Two Sample t-test
data: df_FND$rrst_sensitivity[df_FND$psychotropic_medication == "1"] and df_FND$rrst_sensitivity[df_FND$psychotropic_medication == "0"]
t = -2.1273, df = 37.077, p-value = 0.04011
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-2.00230723 -0.04882873
sample estimates:
mean of x mean of y
2.309574 3.335142 res<-t.test(df_FND$metascore[df_FND$psychotropic_medication == "1"], df_FND$metascore[df_FND$psychotropic_medication=="0"])
res #
Welch Two Sample t-test
data: df_FND$metascore[df_FND$psychotropic_medication == "1"] and df_FND$metascore[df_FND$psychotropic_medication == "0"]
t = -1.4203, df = 39.126, p-value = 0.1634
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.07204124 0.01260121
sample estimates:
mean of x mean of y
0.5533079 0.5830279 res<-t.test(df_FND$sdq[df_FND$psychotropic_medication == "1"], df_FND$sdq[df_FND$psychotropic_medication=="0"])
res # sig difference across groups with or without intake of psychotropic medication: t = 2.6337, p-value = 0.0129
Welch Two Sample t-test
data: df_FND$sdq[df_FND$psychotropic_medication == "1"] and df_FND$sdq[df_FND$psychotropic_medication == "0"]
t = 2.6337, df = 32.012, p-value = 0.0129
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
2.584876 20.230913
sample estimates:
mean of x mean of y
46.15789 34.75000 # differences of sex in our predictor and outcome variable?
res<-t.test(df_FND$rrst_sensitivity[df_FND$sex=="1"], df_FND$rrst_sensitivity[df_FND$sex=="2"])
res #
Welch Two Sample t-test
data: df_FND$rrst_sensitivity[df_FND$sex == "1"] and df_FND$rrst_sensitivity[df_FND$sex == "2"]
t = -1.1197, df = 15.303, p-value = 0.2801
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-2.0795758 0.6455268
sample estimates:
mean of x mean of y
2.681884 3.398908 res<-t.test(df_FND$metascore[df_FND$sex=="1"], df_FND$metascore[df_FND$sex=="2"])
res #
Welch Two Sample t-test
data: df_FND$metascore[df_FND$sex == "1"] and df_FND$metascore[df_FND$sex == "2"]
t = 1.1436, df = 16.994, p-value = 0.2687
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.02470859 0.08319037
sample estimates:
mean of x mean of y
0.5780561 0.5488152 res<-t.test(df_FND$sdq[df_FND$sex=="1"], df_FND$sdq[df_FND$sex=="2"])
res #
Welch Two Sample t-test
data: df_FND$sdq[df_FND$sex == "1"] and df_FND$sdq[df_FND$sex == "2"]
t = 0.33612, df = 17.95, p-value = 0.7407
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-9.402255 12.982900
sample estimates:
mean of x mean of y
40.29032 38.50000 Seeing that the intake of psychotropic medication is both correlated to our outcome and predictor variables, we have to include this in our linear regression. We add it in an interaction term with group, as it is correlated with group, and then we run it on the full model (including both groups). For the other variables (sex, age, bmi) we do not have to include those, as they are not statistically correlated to both variables. For affective symptoms - there is a trend, hence we run a supplemantary analysis for including these as well.
4.4) Linear Regression
# RRST sensitivity
sdq_lm <- lm(formula = sdq ~ rrst_sensitivity*group + psychotropic_medication,
data = df)
summary(sdq_lm)
Call:
lm(formula = sdq ~ rrst_sensitivity * group + psychotropic_medication,
data = df)
Residuals:
Min 1Q Median 3Q Max
-23.741 -4.157 -1.555 3.758 38.900
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.4719 4.8244 4.658 1.18e-05 ***
rrst_sensitivity 0.3626 1.2835 0.282 0.778274
group1 21.9827 5.8059 3.786 0.000286 ***
psychotropic_medication1 7.2866 2.8532 2.554 0.012460 *
rrst_sensitivity:group1 -3.0980 1.5711 -1.972 0.051916 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.575 on 84 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.4898, Adjusted R-squared: 0.4655
F-statistic: 20.16 on 4 and 84 DF, p-value: 1.146e-11confint(sdq_lm, level = 0.95) 2.5 % 97.5 %
(Intercept) 12.878160 32.06567055
rrst_sensitivity -2.189855 2.91497870
group1 10.436956 33.52837270
psychotropic_medication1 1.612716 12.96052136
rrst_sensitivity:group1 -6.222352 0.02630844 # Save as a text file
lm_summary <- capture.output(summary(sdq_lm))
lm_confint <- capture.output(confint(sdq_lm, level = 0.95))
lm_results <- c("### Linear Model Summary ###", lm_summary,
"", "### Confidence Intervals ###", lm_confint)
writeLines(lm_results, "4.4_lm_rrst_med.txt")
sdq_lm_affective <- lm(formula = sdq ~ rrst_sensitivity*group + psychotropic_medication + anx_dep_SUM,
data = df)
summary(sdq_lm_affective) #
Call:
lm(formula = sdq ~ rrst_sensitivity * group + psychotropic_medication +
anx_dep_SUM, data = df)
Residuals:
Min 1Q Median 3Q Max
-24.061 -4.445 -1.013 3.018 38.181
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.57599 5.06195 3.867 0.000218 ***
rrst_sensitivity -0.03933 1.29067 -0.030 0.975765
group1 19.12577 5.97891 3.199 0.001954 **
psychotropic_medication1 5.04848 3.11010 1.623 0.108327
anx_dep_SUM 0.10520 0.06154 1.710 0.091094 .
rrst_sensitivity:group1 -2.67953 1.57261 -1.704 0.092146 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.468 on 83 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.5071, Adjusted R-squared: 0.4775
F-statistic: 17.08 on 5 and 83 DF, p-value: 1.39e-11confint(sdq_lm_affective, level = 0.95) 2.5 % 97.5 %
(Intercept) 9.50798058 29.6440034
rrst_sensitivity -2.60641766 2.5277619
group1 7.23395740 31.0175789
psychotropic_medication1 -1.13737894 11.2343361
anx_dep_SUM -0.01719769 0.2276056
rrst_sensitivity:group1 -5.80738085 0.4483295 # Save as a text file
lm_summary <- capture.output(summary(sdq_lm_affective))
lm_confint <- capture.output(confint(sdq_lm_affective, level = 0.95))
lm_results <- c("### Linear Model Summary ###", lm_summary,
"", "### Confidence Intervals ###", lm_confint)
writeLines(lm_results, "4.4_lm_rrst_med_aff.txt")
# metascore
sdq_lm <- lm(formula = sdq ~ metascore*group + psychotropic_medication,
data = df)
summary(sdq_lm) #
Call:
lm(formula = sdq ~ metascore * group + psychotropic_medication,
data = df)
Residuals:
Min 1Q Median 3Q Max
-25.427 -3.991 -1.279 3.525 38.022
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 18.253 21.722 0.840 0.4031
metascore 9.522 37.609 0.253 0.8007
group1 54.342 25.040 2.170 0.0328 *
psychotropic_medication1 7.997 2.760 2.897 0.0048 **
metascore:group1 -73.283 43.279 -1.693 0.0941 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.56 on 84 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.4914, Adjusted R-squared: 0.4672
F-statistic: 20.29 on 4 and 84 DF, p-value: 1.005e-11confint(sdq_lm, level = 0.95) 2.5 % 97.5 %
(Intercept) -24.944466 61.44956
metascore -65.266950 84.31116
group1 4.546173 104.13700
psychotropic_medication1 2.508047 13.48622
metascore:group1 -159.348639 12.78268 # Save as a text file
lm_summary <- capture.output(summary(sdq_lm))
lm_confint <- capture.output(confint(sdq_lm, level = 0.95))
lm_results <- c("### Linear Model Summary ###", lm_summary,
"", "### Confidence Intervals ###", lm_confint)
writeLines(lm_results, "4.4lm_meta_med.txt")
sdq_lm_affective <- lm(formula = sdq ~ metascore*group + psychotropic_medication + anx_dep_SUM,
data = df)
summary(sdq_lm_affective) #
Call:
lm(formula = sdq ~ metascore * group + psychotropic_medication +
anx_dep_SUM, data = df)
Residuals:
Min 1Q Median 3Q Max
-23.280 -4.199 -1.335 3.362 34.590
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.41127 21.50568 0.717 0.4756
metascore 6.75580 37.16375 0.182 0.8562
group1 51.34935 24.77943 2.072 0.0413 *
psychotropic_medication1 5.78845 2.99396 1.933 0.0566 .
anx_dep_SUM 0.10755 0.06037 1.782 0.0785 .
metascore:group1 -70.57624 42.75691 -1.651 0.1026
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.439 on 83 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.5102, Adjusted R-squared: 0.4806
F-statistic: 17.29 on 5 and 83 DF, p-value: 1.088e-11confint(sdq_lm_affective, level = 0.95) 2.5 % 97.5 %
(Intercept) -27.36266743 58.1852148
metascore -67.16139890 80.6729973
group1 2.06405554 100.6346377
psychotropic_medication1 -0.16641925 11.7433184
anx_dep_SUM -0.01252055 0.2276273
metascore:group1 -155.61801619 14.4655269 # Save as a text file
lm_summary <- capture.output(summary(sdq_lm_affective))
lm_confint <- capture.output(confint(sdq_lm_affective, level = 0.95))
lm_results <- c("### Linear Model Summary ###", lm_summary,
"", "### Confidence Intervals ###", lm_confint)
writeLines(lm_results, "4.4_lm_meta_med_aff.txt")4.5) SUPPLEMENTARY Regression - separate per group and controlling for ANX-DEP sum scores
# RRST sensitivity FND
sdq_lm_FND <- lm(formula = sdq ~ rrst_sensitivity + psychotropic_medication,
data = df_FND)
summary(sdq_lm_FND)
Call:
lm(formula = sdq ~ rrst_sensitivity + psychotropic_medication,
data = df_FND)
Residuals:
Min 1Q Median 3Q Max
-24.141 -8.096 -2.786 6.417 37.778
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 43.391 5.165 8.401 2.28e-10 ***
rrst_sensitivity -2.591 1.323 -1.959 0.0571 .
psychotropic_medication1 8.751 4.263 2.053 0.0467 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.16 on 40 degrees of freedom
Multiple R-squared: 0.2279, Adjusted R-squared: 0.1893
F-statistic: 5.903 on 2 and 40 DF, p-value: 0.005668confint(sdq_lm_FND, level = 0.95) 2.5 % 97.5 %
(Intercept) 32.9527048 53.82919851
rrst_sensitivity -5.2640573 0.08229745
psychotropic_medication1 0.1358901 17.36565231sdq_lm_affective_FND <- lm(formula = sdq ~ rrst_sensitivity + psychotropic_medication + anx_dep_SUM,
data = df_FND)
summary(sdq_lm_affective_FND) #
Call:
lm(formula = sdq ~ rrst_sensitivity + psychotropic_medication +
anx_dep_SUM, data = df_FND)
Residuals:
Min 1Q Median 3Q Max
-24.312 -8.076 -2.783 5.705 37.430
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 39.34986 7.81643 5.034 1.12e-05 ***
rrst_sensitivity -2.60166 1.33144 -1.954 0.0579 .
psychotropic_medication1 6.88075 5.06997 1.357 0.1825
anx_dep_SUM 0.07699 0.11120 0.692 0.4928
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.25 on 39 degrees of freedom
Multiple R-squared: 0.2373, Adjusted R-squared: 0.1786
F-statistic: 4.044 on 3 and 39 DF, p-value: 0.0135confint(sdq_lm_affective_FND, level = 0.95) 2.5 % 97.5 %
(Intercept) 23.5396403 55.16007691
rrst_sensitivity -5.2947400 0.09142542
psychotropic_medication1 -3.3742320 17.13572965
anx_dep_SUM -0.1479358 0.30190731# RRST sensitivity HC
sdq_lm_HC <- lm(formula = sdq ~ rrst_sensitivity + psychotropic_medication,
data = df_HC)
summary(sdq_lm_HC)
Call:
lm(formula = sdq ~ rrst_sensitivity + psychotropic_medication,
data = df_HC)
Residuals:
Min 1Q Median 3Q Max
-4.328 -3.163 -1.072 1.924 11.822
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 24.45286 2.05771 11.884 3.57e-15 ***
rrst_sensitivity -0.09285 0.54018 -0.172 0.864
psychotropic_medication1 -0.76881 2.93689 -0.262 0.795
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.865 on 43 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.001805, Adjusted R-squared: -0.04462
F-statistic: 0.03888 on 2 and 43 DF, p-value: 0.9619confint(sdq_lm_HC, level = 0.95) 2.5 % 97.5 %
(Intercept) 20.303100 28.6026252
rrst_sensitivity -1.182222 0.9965168
psychotropic_medication1 -6.691613 5.1539933sdq_lm_affective_HC <- lm(formula = sdq ~ rrst_sensitivity + psychotropic_medication + anx_dep_SUM,
data = df_HC)
summary(sdq_lm_affective_HC) #
Call:
lm(formula = sdq ~ rrst_sensitivity + psychotropic_medication +
anx_dep_SUM, data = df_HC)
Residuals:
Min 1Q Median 3Q Max
-6.1591 -2.8961 -0.4383 1.9792 8.4374
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 20.46811 2.17139 9.426 6.35e-12 ***
rrst_sensitivity -0.45618 0.49392 -0.924 0.3610
psychotropic_medication1 -1.35880 2.62913 -0.517 0.6080
anx_dep_SUM 0.12607 0.03657 3.447 0.0013 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.453 on 42 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.222, Adjusted R-squared: 0.1664
F-statistic: 3.994 on 3 and 42 DF, p-value: 0.0137confint(sdq_lm_affective_HC, level = 0.95) 2.5 % 97.5 %
(Intercept) 16.08607859 24.8501447
rrst_sensitivity -1.45294757 0.5405941
psychotropic_medication1 -6.66459477 3.9470002
anx_dep_SUM 0.05226925 0.1998698# metascore FND
sdq_lm_FND <- lm(formula = sdq ~ metascore + psychotropic_medication,
data = df_FND)
summary(sdq_lm_FND) #
Call:
lm(formula = sdq ~ metascore + psychotropic_medication, data = df_FND)
Residuals:
Min 1Q Median 3Q Max
-25.720 -7.438 -3.396 8.541 37.124
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 70.469 17.675 3.987 0.000277 ***
metascore -61.264 29.966 -2.044 0.047528 *
psychotropic_medication1 9.587 4.123 2.326 0.025200 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.11 on 40 degrees of freedom
Multiple R-squared: 0.2339, Adjusted R-squared: 0.1956
F-statistic: 6.106 on 2 and 40 DF, p-value: 0.004851confint(sdq_lm_FND, level = 0.95) 2.5 % 97.5 %
(Intercept) 34.747037 106.190457
metascore -121.827371 -0.701054
psychotropic_medication1 1.255103 17.919140sdq_lm_affective_FND <- lm(formula = sdq ~ metascore + psychotropic_medication + anx_dep_SUM,
data = df_FND)
summary(sdq_lm_affective_FND) #
Call:
lm(formula = sdq ~ metascore + psychotropic_medication + anx_dep_SUM,
data = df_FND)
Residuals:
Min 1Q Median 3Q Max
-24.069 -7.355 -2.813 7.833 34.746
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 66.54227 18.59586 3.578 0.000944 ***
metascore -61.78579 30.15574 -2.049 0.047244 *
psychotropic_medication1 7.64265 4.94639 1.545 0.130400
anx_dep_SUM 0.07988 0.11073 0.721 0.474958
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.19 on 39 degrees of freedom
Multiple R-squared: 0.244, Adjusted R-squared: 0.1858
F-statistic: 4.195 on 3 and 39 DF, p-value: 0.0115confint(sdq_lm_affective_FND, level = 0.95) 2.5 % 97.5 %
(Intercept) 28.9286047 104.1559436
metascore -122.7815257 -0.7900468
psychotropic_medication1 -2.3623702 17.6476663
anx_dep_SUM -0.1440951 0.3038646# metascore HC
sdq_lm_HC <- lm(formula = sdq ~ metascore + psychotropic_medication,
data = df_HC)
summary(sdq_lm_HC) #
Call:
lm(formula = sdq ~ metascore + psychotropic_medication, data = df_HC)
Residuals:
Min 1Q Median 3Q Max
-4.280 -3.087 -1.230 2.018 11.924
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 21.7955 8.8443 2.464 0.0178 *
metascore 4.0164 15.2906 0.263 0.7941
psychotropic_medication1 -0.5258 2.8134 -0.187 0.8526
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.864 on 43 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.002719, Adjusted R-squared: -0.04367
F-statistic: 0.05863 on 2 and 43 DF, p-value: 0.9431confint(sdq_lm_HC, level = 0.95) 2.5 % 97.5 %
(Intercept) 3.959307 39.631774
metascore -26.820060 34.852853
psychotropic_medication1 -6.199508 5.147904sdq_lm_affective_HC <- lm(formula = sdq ~ metascore + psychotropic_medication + anx_dep_SUM,
data = df_HC)
summary(sdq_lm_affective_HC) #
Call:
lm(formula = sdq ~ metascore + psychotropic_medication + anx_dep_SUM,
data = df_HC)
Residuals:
Min 1Q Median 3Q Max
-5.935 -2.825 -0.655 2.306 8.960
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17.65519 8.07978 2.185 0.03451 *
metascore 2.52815 13.80551 0.183 0.85558
psychotropic_medication1 -0.54217 2.53876 -0.214 0.83193
anx_dep_SUM 0.11865 0.03609 3.287 0.00205 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.487 on 42 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.2068, Adjusted R-squared: 0.1501
F-statistic: 3.65 on 3 and 42 DF, p-value: 0.01995confint(sdq_lm_affective_HC, level = 0.95) 2.5 % 97.5 %
(Intercept) 1.3495342 33.9608409
metascore -25.3324965 30.3887882
psychotropic_medication1 -5.6655976 4.5812588
anx_dep_SUM 0.0458067 0.19148525) Modeling Accuracy Interaction
R-script provided from Reviewer 3 and adapted to our data; thus acknolwedgment to them
library(tidyverse)── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ forcats 1.0.0 ✔ tibble 3.2.1
✔ purrr 1.0.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ rstatix::filter() masks dplyr::filter(), stats::filter()
✖ dplyr::lag() masks stats::lag()
✖ car::recode() masks dplyr::recode()
✖ purrr::some() masks car::some()
✖ Hmisc::src() masks dplyr::src()
✖ Hmisc::summarize() masks dplyr::summarize()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(glmmTMB)Warning in check_dep_version(dep_pkg = "TMB"): package version mismatch:
glmmTMB was built with TMB package version 1.9.17
Current TMB package version is 1.9.10
Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)library(ggeffects)
library(patchwork)
library(sjPlot)#refugeeswelcomelibrary(ggplot2)
library(ggthemes)
# Example RRST analysis script: modeling confidence data
# Assumes 'rrst_data' contains:
# - confidence: trial-level confidence ratings (scaled 0-1)
# - stimLevel: stimulus intensity (continuous or ordinal)
# - ResponseAccuracy: binary (0 = incorrect, 1 = correct)
# - group: factor (e.g., "FND", "Control")
# - SDQ_score: continuous clinical variable (e.g., dissociation)
# - pcode: participant ID
rrst_confidence_analysis <- readxl::read_xlsx( "/Users/nataschastoffel/Documents/GitHub/interoception_NS/data/processed/05-2025/rrst_confidence_analysis.xlsx", na = c("NaN", "NA", ""))
# data is in wide format, so we have to change to long
library(tidyverse)
rrst_data <- rrst_confidence_analysis %>%
pivot_longer(
cols = matches("^(accuracy|stim|conf)_trial\\d+$"),
names_to = c(".value", "trial"),
names_pattern = "(accuracy|stim|conf)_trial(\\d+)"
) %>%
rename(
ResponseAccuracy = accuracy,
stimLevel = stim,
confidence = conf
) %>%
mutate(
trial = as.integer(trial),
ResponseAccuracy = as.integer(ResponseAccuracy), # binary
stimLevel = as.numeric(stimLevel), # continuous
confidence = as.numeric(confidence), # continuous
group = as.factor(group),
sdq = as.numeric(sdq),
pcode = as.character(pcode) )
# Note: Consider including (1 | stimLevel) as a random effect to control for stimulus-level variance if needed
rrst_data <- rrst_data %>%
mutate(
ResponseAccuracy = factor(ResponseAccuracy, levels = c(1, 0), labels = c("Correct", "Incorrect")),
group = factor(group), # Make sure group is also a factor if needed
pcode = as.character(pcode)
)5.1) Interaction with GROUP
# Model: Confidence ~ Stimulus × Accuracy × Group + (1|pcode)
conf_fit <- glmmTMB(confidence ~ stimLevel * ResponseAccuracy * group + (1 | pcode) + (1 | stimLevel),
data = rrst_data,
family = ordbeta(),
start = list(psi = c(0, 1)))
summary(conf_fit) Family: ordbeta ( logit )
Formula:
confidence ~ stimLevel * ResponseAccuracy * group + (1 | pcode) +
(1 | stimLevel)
Data: rrst_data
AIC BIC logLik -2*log(L) df.resid
3970.6 4053.9 -1972.3 3944.6 4470
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
pcode (Intercept) 0.37135 0.6094
stimLevel (Intercept) 0.04017 0.2004
Number of obs: 4483, groups: pcode, 91; stimLevel, 17
Dispersion parameter for ordbeta family (): 3.19
Conditional model:
Estimate Std. Error z value
(Intercept) 0.22137 0.29777 0.743
stimLevel 0.03737 0.02184 1.711
ResponseAccuracyIncorrect 0.10508 0.28360 0.370
groupHC -0.42289 0.30412 -1.391
stimLevel:ResponseAccuracyIncorrect -0.04768 0.02088 -2.283
stimLevel:groupHC 0.03203 0.02027 1.580
ResponseAccuracyIncorrect:groupHC 0.26384 0.36591 0.721
stimLevel:ResponseAccuracyIncorrect:groupHC -0.01864 0.02769 -0.673
Pr(>|z|)
(Intercept) 0.4572
stimLevel 0.0871 .
ResponseAccuracyIncorrect 0.7110
groupHC 0.1644
stimLevel:ResponseAccuracyIncorrect 0.0224 *
stimLevel:groupHC 0.1141
ResponseAccuracyIncorrect:groupHC 0.4709
stimLevel:ResponseAccuracyIncorrect:groupHC 0.5008
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Model summary table
tab_model(conf_fit,
pred.labels = c("Intercept",
"Stimulus Intensity",
"Response Accuracy",
"Group (FND vs Control)",
"Stim × Accuracy",
"Stim × Group",
"Accuracy × Group",
"Stim × Accuracy × Group"),
dv.labels = "Confidence",
file = "RRST_confidence_model.doc")| Confidence | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| Intercept | 1.25 | 0.70 – 2.24 | 0.457 |
| Stimulus Intensity | 1.04 | 0.99 – 1.08 | 0.087 |
| Response Accuracy | 1.11 | 0.64 – 1.94 | 0.711 |
| Group (FND vs Control) | 0.66 | 0.36 – 1.19 | 0.164 |
| Stim × Accuracy | 0.95 | 0.92 – 0.99 | 0.022 |
| Stim × Group | 1.03 | 0.99 – 1.07 | 0.114 |
| Accuracy × Group | 1.30 | 0.64 – 2.67 | 0.471 |
| Stim × Accuracy × Group | 0.98 | 0.93 – 1.04 | 0.501 |
| Random Effects | |||
| σ2 | 0.14 | ||
| τ00 pcode | 0.37 | ||
| τ00 stimLevel | 0.04 | ||
| ICC | 0.75 | ||
| N pcode | 91 | ||
| N stimLevel | 17 | ||
| Observations | 4483 | ||
| Marginal R2 / Conditional R2 | 0.106 / 0.775 | ||
# Define palette
color_palette <- c("Correct" = "#1b9e77", "Incorrect" = "#d95f02")
# Plot 1: Group × Accuracy Interaction
plot_group_accuracy <- plot_model(conf_fit, type = "pred", terms = c("group", "ResponseAccuracy")) +
scale_color_manual(values = color_palette) +
labs(title = "A) Group and Accuracy Interaction",
x = "Group",
y = "Predicted Confidence") +
theme_minimal(base_size = 10) +
theme(legend.position = "top",
plot.title = element_text(face = "bold", hjust = 0.5),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8))
# Plot 2: Stimulus × Accuracy × Group Interaction
plot_stim_acc_group <- plot_model(conf_fit, type = "pred", terms = c("stimLevel", "ResponseAccuracy", "group")) +
scale_color_manual(values = color_palette) +
scale_fill_manual(values = color_palette) +
labs(title = "B) Stimulus Intensity × Accuracy × Group Interaction",
x = "Stimulus Intensity",
y = "Predicted Confidence") +
theme_minimal(base_size = 10) +
theme(legend.position = "top",
plot.title = element_text(face = "bold", hjust = 0.5),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8))
# Combine plots
combined_plot <- plot_group_accuracy / plot_stim_acc_group +
plot_layout(guides = 'collect') & theme(legend.position = "bottom")
print(combined_plot)
ggsave("S2_Accuracy_interactions_GROUP.tiff",
plot = combined_plot,
dpi = 600,
width = 6,
height = 4,
units = "in",
device = "tiff",
compression = "lzw")5.2) Interaction with SDQ
# Model: Confidence ~ Stimulus × Accuracy × Group + (1|pcode)
conf_fit <- glmmTMB(confidence ~ stimLevel * ResponseAccuracy * sdq + (1 | pcode) + (1 | stimLevel),
data = rrst_data,
family = ordbeta(),
start = list(psi = c(0, 1)))
summary(conf_fit) Family: ordbeta ( logit )
Formula:
confidence ~ stimLevel * ResponseAccuracy * sdq + (1 | pcode) +
(1 | stimLevel)
Data: rrst_data
AIC BIC logLik -2*log(L) df.resid
3930.6 4013.9 -1952.3 3904.6 4470
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
pcode (Intercept) 0.36705 0.6058
stimLevel (Intercept) 0.04832 0.2198
Number of obs: 4483, groups: pcode, 91; stimLevel, 17
Dispersion parameter for ordbeta family (): 3.22
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.0526862 0.4103804 -2.565 0.01031
stimLevel 0.1552625 0.0288582 5.380 7.44e-08
ResponseAccuracyIncorrect 1.1060807 0.4739671 2.334 0.01961
sdq 0.0331362 0.0115360 2.872 0.00407
stimLevel:ResponseAccuracyIncorrect -0.1499379 0.0369706 -4.056 5.00e-05
stimLevel:sdq -0.0031433 0.0007546 -4.165 3.11e-05
ResponseAccuracyIncorrect:sdq -0.0232996 0.0142408 -1.636 0.10181
stimLevel:ResponseAccuracyIncorrect:sdq 0.0025795 0.0010877 2.372 0.01772
(Intercept) *
stimLevel ***
ResponseAccuracyIncorrect *
sdq **
stimLevel:ResponseAccuracyIncorrect ***
stimLevel:sdq ***
ResponseAccuracyIncorrect:sdq
stimLevel:ResponseAccuracyIncorrect:sdq *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Model summary table
tab_model(conf_fit,
pred.labels = c("Intercept",
"Stimulus Intensity",
"Response Accuracy",
"SDQ",
"Stim × Accuracy",
"Stim × SDQ",
"Accuracy × SDQ",
"Stim × Accuracy × SDQ"),
dv.labels = "Confidence",
file = "RRST_confidence_model.doc")| Confidence | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| Intercept | 0.35 | 0.16 – 0.78 | 0.010 |
| Stimulus Intensity | 1.17 | 1.10 – 1.24 | <0.001 |
| Response Accuracy | 3.02 | 1.19 – 7.65 | 0.020 |
| SDQ | 1.03 | 1.01 – 1.06 | 0.004 |
| Stim × Accuracy | 0.86 | 0.80 – 0.93 | <0.001 |
| Stim × SDQ | 1.00 | 1.00 – 1.00 | <0.001 |
| Accuracy × SDQ | 0.98 | 0.95 – 1.00 | 0.102 |
| Stim × Accuracy × SDQ | 1.00 | 1.00 – 1.00 | 0.018 |
| Random Effects | |||
| σ2 | 0.14 | ||
| τ00 pcode | 0.37 | ||
| τ00 stimLevel | 0.05 | ||
| ICC | 0.75 | ||
| N pcode | 91 | ||
| N stimLevel | 17 | ||
| Observations | 4483 | ||
| Marginal R2 / Conditional R2 | 0.129 / 0.784 | ||
# Define palette
color_palette <- c("Correct" = "#1b9e77", "Incorrect" = "#d95f02")
# Get predicted values so that the confidence inteval sare correcelty colored
library(ggeffects)
pred_sdq_accuracy <- ggpredict(conf_fit, terms = c("sdq", "ResponseAccuracy"))
# Plot 1: SDQ × Accuracy with matched line and ribbon colors
plot_sdq_accuracy <- ggplot(pred_sdq_accuracy, aes(x = x, y = predicted, color = group, fill = group)) +
geom_line(size = 1.2) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2, color = NA) +
scale_color_manual(values = color_palette, name = "ResponseAccuracy") +
scale_fill_manual(values = color_palette, name = "ResponseAccuracy") +
labs(title = "A) SDQ and Accuracy Interaction",
x = "SDQ",
y = "Predicted Confidence") +
theme_minimal(base_size = 10) +
theme(legend.position = "top",
plot.title = element_text(face = "bold", hjust = 0.5),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8))
# Plot 2: Stimulus × Accuracy × SDQ Interaction
plot_stim_acc_sdq <- plot_model(conf_fit, type = "pred", terms = c("stimLevel", "ResponseAccuracy", "sdq")) +
scale_color_manual(values = color_palette) +
scale_fill_manual(values = color_palette) +
labs(title = "B) Stimulus Intensity × Accuracy × SDQ Interaction",
x = "Stimulus Intensity",
y = "Predicted Confidence") +
theme_minimal(base_size = 10) +
theme(legend.position = "top",
plot.title = element_text(face = "bold", hjust = 0.5),
axis.title = element_text(size = 10),
axis.text = element_text(size = 8))
# Combine plots
combined_plot <- plot_sdq_accuracy / plot_stim_acc_sdq +
plot_layout(guides = 'collect') & theme(legend.position = "bottom")
print(combined_plot)
ggsave("S3_Accuracy_interactions_SDQ.tiff",
plot = combined_plot,
dpi = 600,
width = 6,
height = 4,
units = "in",
device = "tiff",
compression = "lzw")